Non-Linear Elasticity in Aluminium Alloys 
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SUMMARY: Tests have been carried out on round bar specimens of various 
strong aluminium alloys (unclad) in the fully heat-treated condition in both 
tensile and compressive loading. A careful study of the results establishes 
beyond reasonable doubt that there is a progressive change in the tangent 
elastic modulus with stress. The variation is appreciably greater than that 
predicted by the consideration of second-order effects in the conventional stress 
and strain representation. It is found that the modulus decreases with tension 
and increases with compression by about 5 per cent at 30 tons/in.?, the form 
of the variation being a continuous function. Tests on two alloys of aluminium 
other than those of the normal zinc-bearing and the copper-bearing types were 
also carried out. The first of these was a simple fully heat-treated binary 
alloy of 4 per cent copper which was very similar in behaviour to D.T.D.364; 
the second was a 5 per cent alloy (N.S.6-4H) which, after its initial loading, 
showed only very slight elastic modulus variation in the subsequent loading 
and unloading lines. An additional test to compare naturally and artificially 
aged D.T.D.364 did not reveal any difference in the degree of modulus variation. 


In view of the present experimental findings, the use of the “second” 
modulus line in the proof stress determination of unclad strong alloys should 
be adopted generally, if practicable. A high level of internal stress can cause 
slight errors in this method because the elastic limit is then prematurely reached. 
As an occasional check on the method when applied to routine mechanical 
testing, the elastic component of strain can be found for the last point on the 
loading line by returning to the zeroing load. 


In the Appendix the difference between true and nominal stress-strain 
relationships is calculated. The meaning and theory of modulus variation is 
also discussed and aluminium alloys are compared with a high-tensile steel. 


A simple expression for the form of the variation is proposed. 


1. Introduction 


The non-linear elasticity of metallic materials is a direct contradiction of the 
simple laws which are commonly regarded to hold true. In attempting to decide 
for any material whether or not there is non-linear elastic behaviour, it is important 
to eliminate a number of factors which can, and do, modify the shape of a stress- 
strain curve such as is obtained in a proof stress determination. If a particular 
material does exhibit a genuine curvature of the purely elastic portion of its stress- 
strain curve when all these factors have been eliminated, then the result has 
far-reaching implications in several different fields. 


*The work described in this paper was done while the author was a research engineer with The 
British Aluminium Company. 
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Firstly, those concerned with routine mechanical testing must acknowledge 
that a curved elastic line can exist and must adapt their procedure accordingly. At 
present, it is fairly generally recognised that unclad strong aluminium alloys show 
what is commonly called the “double modulus” effect. A.I.D. Proof Stress 
Determination Procedure has permitted the use, for D.T.D. 687, of what appears in 
routine testing to be a second straight line on the elastic portion of the stress-strain 
curve, and hence the use of the term “double modulus.” The difference in the 
value of proof stress found by using the “ first” and “ second ” lines is between two 
and four per cent. If it can be shown that both these lines are merely tangents to a 
continuous elastic curve, then the “second” line should be used exclusively for 
unclad materials, and in many cases the difference between proof stress and 
permanent set stress would then disappear. 


Secondly, designers could use accurate tangent modulus figures pertaining to 
the stress at which they seek to operate. Thirdly, metallurgists should direct their 
attention to determining the factors which influence the behaviour. 


It is important at this point to make a clear distinction between the tangent 
modulus of the “stress versus total strain” curve and the corresponding tangent 
modulus of the “ stress versus elastic strain component” curve. Up to the elastic 
limit these quantities are identical, but beyond it the former falls rapidly while the 
latter shows no such sudden change. In the compressive loading of a strut it is the 
total tangent modulus which governs buckling. 


Although the ratio of the amount of compression testing to tension testing is 
very low, it is nevertheless the compression properties of a material which are often 
of prime importance to the aircraft designer. Tests have therefore been carried 
out in the present work to examine the variation, if any, in tangent modulus in 
compressive loading. It is fully realised that axial loading is far more difficult to 
achieve in compression than in tension and that the results of compression tests are 
possibly not so reliable as those of tension tests. However, the form of the modulus 
characteristics for the full stress range gives confidence to the acceptance of the 
compression results. 


NOTATION 
Pp nominal stress 
E,E’ Young’s modulus* and tangent modulus respectively 
E, modulus at zero stress 
1/m  Poisson’s ratio 
e,é conventional strain and natural strain respectively 
A, _ initial cross-sectional area 


A. cross-sectional area when the axial strain is e 


*j.e. the ratio of nominal stress to total elastic strain. 
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2. Preliminary Investigations 
2.1. THE EFFECT OF RESIDUAL STRESS 


A series of proof stress tests was carried out on D.T.D. 687 (Al Zn Mg Cu alloy) 
fully heat-treated plate material in the unclad condition, in connection with an 
investigation into residual stresses, and the extension was measured by a dial gauge 
type of extensometer. Inspection of these curves showed a slight but definite 
curvature within the elastic portion, which indicated an apparent drop in tangent 
modulus with applied stress. Those samples which were known to possess a high 
level of residual stress showed a greater change in slope than did those which were 
virtually stress free. However, in all cases the elastic lines obtained on re-loading 
after a small permanent extension had been produced were the same. These showed 
a change in slope which was the same as those of the unstressed samples on initial 
loading. 


The effect of residual stress on the shape and size of the elasto-plastic fillet, as 
the “ knee ” of the stress-strain curve has been called, is shown in Fig. 1. The core 
of quenched material is normally in tension, and it is this region that reaches its 
elastic limit first. The specimen has, therefore, an appreciably lower elastic limit 


(a) UNSTRESSE 


b) HIGH INTERNAL 
STRESS 


ELASTIC 
LIMIT 
(a) 
Fig. 1. 
Diagram of stress-strain curves for 
D.T.D. 687 in (a) unstressed and 
ELASTIC , (b) high internal stress conditions. 
LIMIT 
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than its unstressed counterpart. For the particular material in the tests, it was found 
that the highly stressed samples had an elastic limit of about 20 tons/in.”, whereas 
those with a very low stress level had an elastic limit of nearly 30 tons/in.’, i.e. very 
near the proof stress. These particular values would seem to indicate that the maxi- 
mum tensile residual stress was about 10 tons/in.? This value agrees well with the 
maximum stress which could be carried elastically immediately after quenching, and 
it is reasonable to assume that this value could persist through precipitation heat 
treatment without a major drop. 


Before leaving residual stress and its effect on the variation in tangent modulus, 
it is interesting to note that, as the size and shape of the elasto-plastic fillet is 
governed by the residual stress level on a macroscopic scale, it is reasonable to 
assume that the very presence of the fillet is due to residual stress. This can be 
either on a limited macroscopic scale, or on a purely microscopic scale, as the 
difference in stresses carried by adjacent grains. 


2.2. THE EFFECT OF CLADDING 


Although the inference of the tests of Section 2.1 is that in the absence of clad- 
ding there is curvature of the loading line, it is interesting to compare the shapes 
of the loading lines for two identical samples from a clad sheet, one with the 
cladding subsequently removed by milling. The test on the clad specimen (D.T.D. 
687) was carried out at the A.I.D. Laboratories, Harefield, as part of routine testing. 
At the author’s request, an identical specimen from the same sheet (0-160 in.) was 
tested with the Al 1% Zn cladding removed, in the same testing machine 
and using the same extensometer. The resulting curves are shown in Fig. 2. A 
modulus line of 10-3 x 10° lb./in.? has been included on the graph, this being the 
common value of tangent modulus in the 3-5 tons/in.? region. Both curves show the 
“double modulus” effect, although, as would be expected, the unclad material is less 
subject to a drop in modulus with applied stress. 


3. Main Tests 


In view of the results of Section 2, further tests to investigate fully the non-linear 
elastic behaviour Of several high-strength alloys were made with the co-operation of 
A.I.D. laboratories. A newly installed 30 tons Amsler machine was used for the 
tensile tests; because of the height of the compression platens of this machine, the 
compression test results were obtained on a 10 tons Avery machine, and one of 
the tensile tests was repeated on the Avery for direct comparison with the Amsler. 
The same 2 in. gauge length Lamb roller extensometer was used for all tests, both 
in tension and compression. 


4. Apparatus 
4.1. TESTING MACHINES 


Both testing machines had been calibrated a short time before the tests. They 
complied with the requirements of B.S.S. 1610 for Grade A testing machines; any 
errors from nominal were insignificant. 
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The machines were used entirely on their 10-ton ranges, except for the one test 
on a sheet specimen for which the 3-ton range of the Amsler was used. 


4.2. EXTENSOMETER 


The 2 in. gauge length Lamb roller extensometer used for the present work was 
the subject of an A.I.D. Report, No. A.I.D./Metro/2 and Supplement. The 
instrument had been carefully and accurately calibrated with slip gauges. It was 
found that a modification to the instrument was necessary, because preliminary work 
showed that there could be loss of motion between the specimen and the plates due 
to the resistance offered by the main clamping springs. This possible source of 
error at large strains was overcome by inserting two pairs of small plates between 
the springs and the main part of the instrument; one of each pair of plates carried 
four steel balls trapped in place by a brass retaining cover, and the other was a plain 
steel plate with a lapped surface. 
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When the modification had been made the extensometer gave a repeatable 
response to within 0-00001 in. over a range of 0-015 in. extension. Using a 150 cm. 
curved scale and telescope in preference to a light source, the mean accuracy of the 
instrument was found to be +0-000000 in. and —0-00001 in. The errors were not 
cumulative over the full range and, considering that readings were rounded off to the 
nearest 0-000005 in., the accuracy and repeatability of the instrument may be said 
to be of a very high order. 


One important factor in the use of an accurate extensometer is the effect of 
thermal expansion on the specimen and extensometer frame. In this connection, the 
close proximity of an investigator to the specimen under test could produce a tem- 
perature rise sufficient to alter appreciably the extensometer reading. The present 
work was carried out in a laboratory free from draughts and in which the ambient 
temperature remained at a sensibly constant value. 


4.3. TEST SPECIMENS 


In all but one case the specimens were round bars with a parallel portion of at 
least 3 in. accurately ground to 0:564+0-0002 in. diameter. The specimen ends 
were threaded with } in. diameter Whitworth thread; all operations were carried 
out between centres. 


The materials fulfilled D.T.D. Specifications and, to avoid destroying or damag- 
ing the specimens after testing, spectrographic analysis only was carried out to check 
the approximate levels of the alloying constituents. The following results were 
obtained : — 


Specification Constituents 
D.T.D. 687 (unclad) Zn approx. 6%; Mg approx. 2°5% 
Cr 0:1%; Cu and Mn (not determined) 
D.T.D. 683 Zn approx. 6%; Mg approx. 2°5% 
Cr 0:1%; Cu and Mn (not determined) 
D.T.D. 363 (i) and (ii) Zn 6-64%; Mg approx. 3%; Cr 0°1%; 
Cu and Mn (not determined) 
D.T.D. 364 (i) and (ii) Cu 3-4%; Si, Fe, Mn, and Mg approx. 
0:5% each. 


In addition to these specimens, one round bar specimen of an alloy of 4% Cu 
with 99:99% Al was tested, together with a sheet specimen of N.S. 6 (5% Mg) alloy 
in the cold-rolled condition. All the heat-treatable materials including the 4% Cu 
were in the fully heat-treated condition, except for the one additional test on 
D.T.D. 364. 


5. Method of Testing 


Before the start of a test, the testing machine was run for a period under 
load to allow the oil to circulate and become warm. The specimen was fitted into 
the jaws of the machine and the Lamb extensometer attached to it. Careful align- 
ment and adjustment of rollers and mirrors was then carried out so that the scale, 
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which was accurately positioned, could be read from a telescope with fine cross 
lines. Zeroing was then carried out at a load corresponding to 1 ton/in.* 


Loading was by increments of 0-25 ton, i.e. 1 ton/in.*, and the scale was read 
to the nearest } mm., i.e. giving a sensitivity of 0-0000025 in. extension. The 
operator responsible for loading allowed a very gradual rise through the required 
value, which, for the round bars, was always a multiple of 0:25 ton. A verbal 
indication was given to the observer at the telescope at the precise moment when 
the load was reached. This method was found to give repeatable results. 


The third member of the team recorded the load and extension values, and 
plotted the strain increments for each unit of stress. When the strain increments 
showed signs of a progressive increase, or alternatively at about 20 tons/in.” stress, 
the specimen was unloaded and readings taken as the load fell through the required 
values. The residual extension or compression at the zeroing load was recorded, 
both as a descending value and, after a short pause, as an ascending one. Re-loading 
values were then taken up to and beyond the previous maximum stress level; a 
second unloading and re-loading usually took place when the permanent strain was 
about 0-05 per cent and again at the proof strain of 0-1 per cent. At this point the 
specimen was in some cases subjected to a “ scragging ” treatment by repeated load- 
ing up to a stress just below the maximum reached. A careful survey of this new 
elastic line was then made, taking both ascending and descending values. 


For compression testing, it was necessary to construct two steel pads with 
central tapped holes to carry the threaded ends of the specimens. The faces of the 
pads were ground square with the axis of the thread so that the external faces of the 
two pads, when screwed on the specimen, would be as nearly parallel as possible. 
In practice there was about 0-002 in. variation in this dimension. 


Tension and compression tests were in no way different, except that the Avery 
machine was used for the compression tests. A preliminary calculation of the 
buckling load for the specimens in compression gave a value of 120,000 Ib. for the 
fixed-ended condition if the load was initially axial. Loading without the extenso- 
meter confirmed that the specimen was stable at 67,200 Ib./in.’, i.e. 30 tons/in.?, 
provided that the proof stress had not been reached. It was decided to restrict the 
compressive stress to this value, as it was possible that the truly fixed-ended 
condition was not realised; non-axial loading might combine with this and with a 
drop in the slope of the stress-strain curve as the specimen became plastic and so 
cause rapid buckling. 


A further difficulty occurs in the compression of a specimen which has 
previously been strained plastically in tension due to the Bauschinger effect. The 
D.T.D. 687 specimen showed an elastic limit of 12 tons/in.? under these conditions, 
but the rate of yielding was slow and it was possible to reach a stress of 24 tons/in.’ 
with only 0:05 per cent residual compressive strain. However, in view of the 
difficulty experienced, a second, unstrained, D.T.D. 363 specimen was tested in com- 
pression and this showed an elastic limit of about 25 tons/in.?, which was also the 
value obtained in tension. A second specimen was also used for the testing of 
D.T.D. 364 in compression and an elastic limit of about 21 tons/in.? was obtained, 
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higher than that for the first specimen in tension (15 tons/in.”). In view of this 
difference, and the need to compare the tensile modulus variation on the two testing 
machines, the second D.T.D. 364 specimen was loaded in tension on the Avery 
machine immediately after receiving 0-05 per cent permanent compressive strain. 
The Bauschinger effect was then extremely pronounced, causing the elastic limit to 
fall to about 7 tons/in.? 


The D.T.D. 683 specimen was not tested in compression, as its tensile behaviour 
had been very similar to the D.T.D. 687. The 4% Cu and 5% Mg specimens 
were also tested only in tension. 


6, Results 
6.1. METHOD OF ANALYSIS 

The analysis of the values obtained in all tests requires a technique which will 
reveal genuine changes in elastic modulus but eliminate random variations in the 


strain increments associated with machine, extensometer and personal errors. The 
method adopted was to plot the strain increments against stress and to draw through 
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these points a line such that the experimental points lay disposed about so that the 
algebraic sum of the errors was a minimum at all stages. The total of the strain 
increments given by the line would then equal the total strain recorded in the test 
as nearly as possible, as the form of the errors in machine, extensometer and 
operators was in no way cumulative. The difference between the plot of experi- 
mental strain increments and those corrected from the calibration figures for the 
testing machine was not significantly different. In all cases the figures were therefore 
plotted without applying corrections. Figure 3 shows a typical plot of strain 
increments; the values are for the initial loading of the D.T.D. 683 specimen. The 
broken line shows the mean curve for the re-loading of the specimen after a small 
permanent strain; this is not significantly different from the initial loading line, 
except in the region of low stress and beyond the elastic limit. 


6.2. RESULTS FOR TENSION 


Using the analysis of the strain increments and converting the resulting curves 
to those of tangent modulus, the initial loading and the re-loading lines from the 
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Fig. 5. 
Stress plotted against tangent modulus for the loading of three specimens, 
D.T.D. 364, 5% Mg and 4% Cu. 


proof strain were obtained in tension for all materials tested, and are shown in 
Figs. 4 and 5. 


6.3. HYSTERESIS 

It was found that some specimens showed little or no hysteresis, i.e. the 
unloading and re-loading lines after a small permanent strain were the same, both in 
shape and position. Any slight looping consisted of two parallel lines joined by 
cusps at the upper and lower stress limits. The three zinc-bearing alloys were of 
this type. The copper-bearing alloy, D.T.D. 364, was different in two respects: the 
width of the hysteresis loop was appreciably greater (3-5 times) and also the 
unloading line was not parallel with the re-loading line. The 4% Cu alloy showed 
the latter effect even more pronouncedly, and the change in slope of the unloading 
line was in a different sense from the re-loading line. Fig. 6 shows the tangent 
modulus curves for the unloading and re-loading of three alloys, D.T.D. 687, 
D.T.D. 364 and the 4% Cu. 


6.4. TANGENT MODULUS VARIATION OVER FULL STRESS RANGE 

Combining the results of tension and compression for D.T.D. 687, D.T.D. 363 
and D.T.D. 364, the characteristics for the full stress range were obtained for both 
loading and unloading. 
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Fig. 6. 
Comparison between unloading and re-loading lines for three alloys. 


The D.T.D. 687 specimen was loaded in compression after tension on different 
machines, the D.T.D. 364 specimen was loaded in tension after compression on the 
same machine, while two D.T.D. 363 specimens were used in different machines. 
The continuity through zero stress was very good in the case of the two zinc-bearing 
alloys, as shown in Fig. 7. The copper-bearing alloy was again different and a rapid 
change in modulus occurred as the stress passed from tension to compression. 


The interesting point regarding the compression range is that the modulus 
increases. This fact is a reassuring one to strut designers. 


7. Discussion of Results 
7.1. GENERAL 


The results of the tests establish beyond reasonable doubt that the elastic moduli 
of strong aluminium alloys are not simple constants. It is possible that the methods 
used could give rise to errors in the absolute value of the modulus characteristics, 
but not to any substantial amount in the general form of the variation. It is worth 
examining the factors which could give rise to a variation for large strains. There 
are two to consider : (a) the difference between true and nominal stress values and 
the corresponding natural and conventional strain, and (b) a “ second-order ” effect 
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in the force/(elastic displacement) ratio of the ultimate structure of the metal or 
alloy. The second of these factors is really beyond the scope of this paper, but is 
considered briefly in Appendix II, where a comparison is made between the results 
obtained in the present work and those for a high-tensile steel. 


The first factor is concerned with the dimensional changes during elastic 
deformation as affecting the measurement of stress and strain. Both the contraction 
of the stressed area and the increase in gauge length would produce a drop in tangent 
modulus on the conventional tensile stress-strain curve if the true “E” remained 
constant. The equations relating these quantities are derived in the Appendix I, 
where the nominal stress p and the tangent modulus E’ are given by 


where E, is the modulus at zero stress, assumed constant, e is the conventional 
strain, and 1/m is Poisson’s ratio. 
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Comparison between predicted and experimental modulus change in tensile loading for 
three alloys. 


Using these equations, the predicted change in nominal tangent modulus was 
found for 1/m=0-3. This cuvre is shown in Fig. 8 for the tensile range of stress 
used in the tests. On the same figure have been included the experimental curves 
for D.T.D.363, D.T.D.364 and D.T.D.687 converted to a percentage drop in tangent 
modulus. Although the sense of the modulus change is found by experiment to be 
the same as that predicted, the magnitude is far greater. An alternative method of 
measuring tangent modulus which is not affected by dimensional change would 
record the difference between the experimental and predicted curves in Fig. 8. Such 
a method is at present in use at the National Physical Laboratory; a technique is 
being used to measure the velocity of ultrasonic pulses through successive reflections 
along a stressed test piece. Preliminary results on an alloy of the D.T.D. 364 type 
have indicated a fall in modulus with tensile stress of a very similar order to that 
given in Fig. 8. (The author is indebted to Mr. H. L. Cox and Mr. M. F. Markham 
of the National Physical Laboratory for this information). 


May 1957 115 


TT 
N 


A. B. WATTS 


[— LONGITUDINAL 
TRANSVERSE 
4 RMA\ 
TRANSVERSE 
<= 
NORMAL 20 


TENSILE sTRESS ( 


3 


ie) Ol O2 0-3 0-4 05 0-6 07 9-0 9-5 10-0 105 


TENSILE STRAIN (PER CENT) uo 
10°L8./IN: 


Fig. 9. 
Tensile stress-strain curves for 2} in. thick D.T.D. 687 plate tested in longitudinal, transverse 
and normal to rolled surface directions. 


7.2. ANISOTROPY 


The existence of tangent modulus variation with applied stress raises a number 
of questions. One of these is whether all directions within a material are equally 
affected, i.e. whether anisotropy is in any way responsible. This point was settled by 
the analysis of tensile results on D.T.D. 687 (unclad) plate 24 in. thick. The A.I.D. 
Mechanical Test Laboratory kindly allowed the values recorded in the testing of this 
material, using Lamb extensometers, to be analysed. The longitudinal and trans- 
verse specimens had been tested using the same instrument as that for the present 
work, while the normal specimen (short transverse) had necessitated a 1 in. 
instrument. Fig. 9 shows the stress-strain curves for the three principal directions 
and the respective derived modulus curves. The form of the variation is identical 
for the three directions, but the absolute values are slightly different. 


7.3. INITIAL LOADING 


Another question raised by the results is whether the initially high modulus 
found for the first loading of the materials shown in Fig. 4 is a true effect, or merely 
arises from extensometry. Some investigators believe that motion is lost between 
the specimen and extensometer points or knife edges during the first loading. An 
apparently higher modulus would result. It is possible that this occurs, for the 
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Tangent modulus for D.T.D. 364 in both naturally and artificially aged conditions. 


complex stress system existing in the specimen beneath the point of contact is 
immediately modified by the applied stress on loading, and further yielding can 
occur. Another explanation of the initially high modulus is the effect of slight 
bending due to non-axial loading. However, the Lamb extensometer is particularly 
free from this source of error, as it records the mean extension of the two sides of 
the specimen. 


By re-setting the extensometer after a small permanent strain, in a test during 
the present work, it was found that the results agreed with the re-loading line and 
not the initial loading. This suggests that the first line is a true one and not the 
result of lost motion due to the embedding of the extensometer contacts. In this 
context, it is worth recording that stiffened aircraft panels have shown a similar 
high resistance to initial deflection in static loading. (From a private communication 
from the Bristol Aeroplane Company). 


7.4. PARTIAL HEAT-TREATMENT 


One aspect not covered by the main tests was a comparison between the effects 
of full and partial heat-treatment. An additional test was therefore carried out with 
two identical specimens of D.T.D. 364, one naturally aged and the other artificially 
aged. On analysing the results, it was found that there was no difference in the 
degree of modulus variation, although the absolute values were slightly different. 
Figure 10 shows the variation obtained for loading for the two specimens, both 
initially and after a small permanent strain. For the artificially aged specimen the 
two loading lines were slightly displaced. 
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8. Influence of Results on Testing Procedure 


One of the principal objects of the present work was to investigate the use of 
“upper” and “lower” modulus lines in the determination of the proof stress of 
strong alloys. The existence of the curved elastic line apparently in all strong alloys 
gives justification for the use of the “ upper ” line, which approximates very closely 
to the continuation of the curved elastic line past the elastic limit. Figure 11 shows 
the comparison between the use of the upper and lower lines for the determination 
of proof stress for unclad D.T.D. 687. The lower line has a modulus value of 
10-3 x 10° Ib./in.*, the upper, 9-65 x 10° Ib./in.*” The points marked with circles are 
the residual (plastic) strain values after unloading from intermediate stresses up to 
and including the proof stress. They have been plotted as offsets to the left of the 
stress-strain curve, and the remaining strain values are therefore the elastic com- 
ponents corresponding to the particular stress points from which the specimen was 
unloaded. 
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Examination of Fig. 11 shows an appreciable difference between the proof stress 
values corresponding to the two lines. In the author’s opinion, the true proof stress 
is the higher value and is the stress producing a permanent (plastic) strain of 0-1 per 
cent. The apparent proof stress as at present defined is the lower value, which is 
obtained by extending the lower part of the initial stress-strain curve. It corresponds 
to a permanent strain of only 0:07 per cent. To avoid distinguishing between 
“ proof stress” and “ permanent set stress,” it would appear desirable to use the 
“ upper ” modulus line for unclad strong alloys when determining proof stress. As 
a check, the residual strain can be found for the final stress value, before the extenso- 
meter is removed, by unloading to the zeroing load. If this value is then plotted 
as an offset from the stress-strain curve, the “upper” line should pass through 
this point. 


The problem of clad materials is a difficult one. No procedure can be satisfac- 
tory, for a check on permanent strain is prevented by the setting up of residual 
stresses on unloading. The cladding is forced into compression when unloaded 
from a high stress, thus leaving the core material in residual tension. Apparently 
the core has suffered permanent strain, but, in fact, this is not necessarily so. It is 
considered that here proof stress determination must of necessity be by an arbitrary 
method. Two alternatives present themselves: (i) the use of the “ lower ” modulus 
line, i.e. the region of the stress-strain curve where core and cladding are both elastic 
is continued to give the graphical construction for proof stress; (ii) the use of a 
modulus line of a fixed value for each material, drawn tangentially to the 
stress-strain curve. 


9. Conclusions 


(i) The existence has been established of a tangent modulus variation with 
stress for fully heat-treated strong aluminium alloys. 


(ii) The variation is in the form of a decrease in modulus with tensile stress 
and an increase with compressive stress. 


(iii) The zinc-bearing alloys, D.T.D. 687, D.T.D. 683 and D.T.D. 363, show 
a similar form of modulus variation in which loading and unloading lines 
are similar. The copper-bearing alloys, typified by D.T.D. 364, show a 
similar loading characteristic to the zinc-bearing, but their unloading 
behaviour is different. 


(iv) Hysteresis looping is not responsible for modulus variation, although it 
may be associated with it, as for D.T.D. 364. 


(v) As well as strong alloys, the.modulus variation has been found in a 
simple binary heat-treatable alloy containing 4% Cu, and a cold-rolled 
non-heat-treatable alloy N.S. 6 - 4H, containing 5% Mg. The former of 
these two resembled D.T.D. 364 in behaviour, while the latter, although 
showing appreciable variation in initial loading, maintained an almost 
constant value on re-loading after a small plastic strain. 
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(vi) The form of the tangent modulus variation was the same in the three 
principal directions of a thick, rolled slab of D.T.D. 687, although the 
absolute values of the moduli were different. 


(vii) An additional test showed that the degree of modulus variation was the 
same for D.T.D. 364 in the naturally aged condition as in the artificially 
aged condition. 


(viii) It is reasonable to use the “ upper” modulus line in the proof stress 
determination of unclad strong alloys. This should be supplemented 
by the unloading of certain specimens after the final extension value has 
been recorded, to determine the elastic and plastic components of the 
corresponding total strain. The “upper” modulus line should pass 
through the point representing the elastic component of strain plotted at 
the appropriate stress value. 


ACKNOWLEDGMENTS 


The author wishes to acknowledge with thanks the assistance given by Aero- 
nautical Inspection Directorate Laboratories, Harefield, in the experimental work, 
and also to thank The British Aluminium Company for permission to publish 


this paper. 
Appendix I 


CALCULATION OF CONVENTIONAL STRESS-STRAIN CURVE AND TANGENT MODULUS 
FOR CONSTANT TRUE “E” 


It is assumed that the transverse strains bear a simple ratio to the axial strain, given 
by Poisson’s ratio, 1/m. Also that the true stress and natural strain increments are 
always related by the same value of E. 


(i) If A, is the initial cross-sectional area, and A, the area when the axial strain 


is e, then 


A, 
=(1 —e/m) 


0 
Lf dL is a small load increment, then the true stress increment is dL/A,, which 
may be expressed as 
aL 
A, 
(ii) The natural strain increment d¢ is related to the conventional strain increment 
by 


(3) 


de 


Now if the ratio of the true stress increment to the natural strain increment is E,, and 
this remains constant, then 


(4) 


dL 
Ade’ 
Using equations (3) and (4), E,= A. 
i.e. =E, (1 - e/m)? (1 +e)-'de. ; (5) 
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For stress-strain curves, the nominal stress p is plotted, and dp=dL/A,. 


The slope of the conventional stress-strain curve is E’, where 


,_ 4p 
E de 
Therefore . . . 


Equations (6) and (7) gives values of p and E’ corresponding exactly to those used 
in the figures showing the experimental results of the tests. 


Appendix II 
THE THEORY OF ELASTIC MODULUS VARIATION 


The work described in this paper was carried out specifically to investigate the 
apparent variation in elastic modulus with stress of high-strength aluminium alloys. 
It has been shown that the variation is appreciably greater than that which could occur 
as the result of dimensional changes of the specimen during elastic straining. It must 
be concluded, therefore, that the modulus variation is a fundamental “ second-order ” 
effect which can well be ignored if the elastic strain is small, but cannot reasonably 
be if the material in question is one of high strength and/or low modulus. 


All physical properties of materials such as coefficients of thermal expansion, 
thermal and electrical conductivity are not constant in value but vary with the level of 
the property in question. It is not unreasonable that the elastic modulus is similarly 
affected. The ionic theory of metals postulates that the metal ions are held in position 
in the lattice by attractive and repulsive forces which are in equilibrium. Force is 
necessary to displace cach ion relative to the next and it is theoretically possible to 
predict the elastic constants for a particular metal. The relationship between force 
and strain will depend on the shape of the curve of energy versus the distance between 
the ions. 


It is obvious that all metals and alloys should show modulus variation to some 
extent, if indeed it is a fundamental property of their ultimate structure. It was 
therefore decided to examine a high-tensile steel which was known to possess a proof 
strain of a similar order to the aluminium alloys which had been examined. It is 
recognised that high-tensile steels do not show any sharp yield as do annealed mild 
steel. A progressive divergence from linearity occurs in the stress-strain curve as loading 
takes place. It is commonly thought that the metal has passed its elastic limit when 
this change in slope becomes apparent. The test carried out in the present work was 
aimed at examining the elastic line of such a steel (EN.25X) when a small permanent 
strain had been caused. 


The specimen was loaded in tension in the same testing machine as the aluminium 
alloy specimens, and the same extensometer was used. The load was raised from a low 
value to near the proof load and back several times after a normal proof stress deter- 
mination had been made. It was hoped that this would prevent creep at the higher 
stresses. A careful survey was then made of the elastic line between 3 and 67 tons/in.’, 
both in ascending and descending values of load. Analysis of the strain increments 
showed a drop in modulus during loading of about 6 per cent (from 29-5 x 10° to 
27-7 x 10° Ib./in.?). There was a small amount of hysteresis on unloading and the 
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Fig. 12. 
Modulus variation plotted against strain for four aluminium 
alloys and a high-tensile steel, EN.25X. 


modulus variation was less; however, on dropping the load below the zeroing valuc 
and increasing to it, there was found to be no permanent strain whatsoever. 


i view of this finding on a high-tensile steel it was decided to replot the curves 
iit Fig. to i: base of strain and to include the steel curve on the same base (Fig. 12). 
The vai ‘ti \ for the 4 per cent CuAl alloy was also plotted on the figure. 


Ther. .. a marked similarity between the curves for the various materials, which 
suggests that a simple mathematical representation can be made of the behaviour. 
Accordingly, an expression of the form 


E’=E,{1- sinh (ze)} 


was tried. It was found that a common value of a gave good agreement for the 
aluminium alloys if the value of 8 was chosen correctly. If e is a percentage strain, 
then a is 1-75 and the corresponding values of 8 are 


Material D.T.D.363 D.T.D.364 D.T.D.687 4% CuAl 
B 0-032 0-041 0-046 0-100 


Similarly, the values for the function to give agreement with the steel variation 
are 2=3-0, B=0-02. 
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On the Measurement of Supersonic Aerofoil 
Drag by Pressure Traverse 
R. E. MEYER 
(Department of Aeronautics, ctiadiateasd of Sydney)* 


SuMMARY: The drag in two-dimensional flow is expressed in terms of the 

distributions of static and stagnation pressures along a traverse line down- 

stream of the aerofoil. Sources of error are discussed with regard to their 

effect on the accuracy of drag measurement in small tunnels at medium 
Mach numbers. 


1. Introduction and Conclusions 
1.1. INTRODUCTION 


The wake traverse is a standard method of measuring aerofoil drag at subsonic 
speeds, but it seems to have been scarcely considered for supersonic flow. At the 
same time, alternative methods often present considerable difficulties; for instance, 
the installation of a strain-gauge balance inside the model seems prohibitively diffi- 
cult in small tunnels. In supersonic flow, moreover, the pressure distributions 
obtained by the traverse afford a possibility to distinguish to some degree the relative 
contributions to the total drag made by different physical effects. 


The purpose of this paper is to attempt an estimate of the main difficulties likely 
to be encountered with the traverse method in two-dimensional flow at medium 
Mach numbers. The work proceeds in two stages, the first of which is the derivation 
of expressions giving the drag in terms of local static pressure, stagnation pressure 
and stagnation temperature on the traverse line. The results (Sections 3 and 5) 
differ significantly from those of Refs. 1 (p. 567) and 2. The second stage is a survey 
of some probable sources of measuring errors and their effect on the value obtained 
for the drag. 


It is difficult to attempt an exhaustive discussion of sources of experimental 
error in advance of actual experiments, and no such attempt will be made here. The 
present results are therefore hardly sufficient for a proper assessment of the potential 
usefulness of the method. At the same time, any serious preliminary discussion of 
experimental difficulties is liable to have a deterrent effect which is misleading. The 
present results suggest a need for some further preliminary work and for a consider- 
able amount of work on the design of experimental details. But they do not 
confirm the pessimism expressed in Ref. 1 (p. 568), and it is felt that they. go far 
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enough to give some substance to the conjecture that no fundamental difficulties 
impede the measurement of drag by traverse in small tunnels with a useful degree 
of accuracy. 


1.2. WAVE DRAG AND WAKE DRAG 


If a rectangular control box with all sides either parallel or normal to the 
undisturbed stream and all corners upstream of the nose shocks (Fig. 1) is chosen, 
then, by the momentum principle, the drag per unit span is 


where a suffix | refers to the undisturbed stream and the integration is only over the 
downstream side of the box, from nose shock to nose shock. The traverse therefore 
extends across regions with flow patterns of quite different type and it is convenient, 
in particular, to distinguish the wake, defined as the region where the stagnation 
pressure differs appreciably from its undisturbed. value, from the rest of the flow, 
where it differs but little. Thus (1) will be written in the form 


Cp=2D/(p,U?c)= ewdy/e+ | crdy/c, ' (2) 


w 


where cw=cy, but the second integral is taken only over the traverse of the wake 
and the first integral, only over the rest of the traverse. For convenience, cw will 
be called the wave drag contribution and cy the friction drag contribution. 


This distinction does not, of course, correspond to a true physical distinction, 
but it comes near enough to it for a rough estimate of the relative orders of magni- 
tude of the integrals in (2) to be obtainable from simple theories, at least for thin 
aerofoils. Such an estimate is useful for a discussion of the influence, on the accuracy 
of cp, of measuring errors affecting, in the first place, only cw or cp. Since the 
flat-plate aerofoil without incidence has no wave drag, it may be expected that for 
any given aerofoil shape there is a thickness ratio, depending on incidence, Mach 
number and Reynolds number, for which the two integrals are equal. The linearised 
theory of inviscid flow and boundary layer theory and experiments for flat plates 
indicate that this occurs, at zero incidence, for thickness ratios among the smallest of 
practical interest. With increasing thickness ratio, and especially with incidence, the 
wave drag rises rapidly. An approximate evaluation of cp from the Pitot traverse 
of the wake of a 35 per cent thick wedge at zero incidence and M=2, reported in 
Ref. 3, suggests that the total friction drag contribution may amount to only between 
10 and 20 per cent of cp in such cases. For relatively thick aerofoils, and for any 
aerofoil at considerable incidence, the accuracy of the wake traverse may therefore 
be expected to be less important than that of the outer traverse; on the other hand, 
the shock losses, and hence the accurate measurement of the stagnation pressure over 
the outer traverse, then become relatively important (Section 3.2). 
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1.3. CONCLUSIONS 


The contribution to cp from the static pressure traverse from nose shock to tail 
shock is important in any case (Section 3.2), and that from tail shock to wake is also 
significant in cases of extensive interaction between tail shock and boundary layer, 
but a static pressure traverse of the wake is not usually required (Section 5.2). The 
main difficulty appears to arise in the determination of the shock positions, and it is 
suggested (Section 4.1) that this may be overcome by the use of a method ensuring 
mutual compensation of the errors made at the individual shocks. 


A theoretical difficulty is encountered in the determination of the contribution 
to Cw from the Pitot pressure traverse between nose shock and tail shock. The 
second-order theory of the flow suffices neither for approximating cw closely in 
terms of stagnation pressure when the aerofoil is thick or at considerable incidence 
(Section 3.2), nor for calculating the stagnation pressure from the measured Pitot and 
static pressures (Section 4.2), and further work on these points appears desirable. 
Some information, however, can be derived from the static pressure traverse, and 
it is suggested that this may suffice in some cases (Section 4.2). 


The Pitot traverse of the wake, on the other hand, may involve experimental 
difficulties. For instance, even a very small “ displacement effect” on the Pitot 
readings, due to the large traverse gradient of total head, causes a serious error in 
Cr (Section 6.2) and, since little appears known on the flow pattern of a Pitot tube in 
a narrow supersonic wake, further experimental work on this point seems desirable 
before the traverse method can be applied with confidence to thin aerofoils at very 
small incidence. 


A significant change in the drag due to heat transfer must often be anticipated 
but it is compensated to a considerable degree (except at relatively low Reynolds 
numbers) by the error incurred in the traverse method by neglecting the variation of 
stagnation temperaure (Section 6.1). Errors due to tunnel interference are discussed 
briefly in Section 2. The effect produced by the minor errors in the measurement of 
the undisturbed Mach number, static pressure and stagnation pressure on the 
accuracy of the total wave and friction drag contributions is easily deduced from the 
expressions for Cw (equation (7)) and cy (equation (14)). 


NOTATION 
c chord 

Cp drag coefficient per unit span 

cr _ friction drag contribution (see equation (2)) 
Cp average of cp across wake 

Cw wave drag contribution (see equation (2)) 

Cy,Cy specific heats 

Cm specific heat per unit mass of model 

D __ distance between the leading edge and the traverse line 
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ratio of Pitot to stagnation pressure (see equation (9)) 
defined by equations (14a) 

defined by equations (19a) 

defined by equation (22) 

Stanton number 

probe length upstream of interaction 

pressure 

measured Pitot pressure 

velocity magnitude 

entropy per unit mass 

time 

velocity components 

terms on the right hand side of equation (7) 

average of w, over traverse segments between nose and tail shocks 
co-ordinates 

average skin-friction coefficient of (one-sided) flat plate 
Mach number of undisturbed stream 

local Mach number 

Reynolds number based on chord and undisturbed stream 
Reynolds number based on the length / 

temperature 

a representative wall temperature 

thermometer temperature 

velocity of undisturbed stream 


incidence 


=(7-D/y 


Colle 

displacement thickness just upstream of the interaction 
relative errors 
stream direction 

defined by equations (14a) 

viscosity 


error magnification factor 


density 
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Pm density of model 
o Prandtl number 


7 thickness ratio 


Suffixes 
0 local stagnation values 


1 values in the undisturbed stream 
m___ refers to the model 
W _ refers to the wake (except in cw) 
w__ refers to the wall 


A dash (’ ) denotes the relative difference between local and undisturbed values, e.g. 
Po =(Po— Por)! Por: 


2. Traverse Position and Tunne] Interference 


The vertical traverse line suggested by the control surface of Fig. 1 implies a 
restriction on the aerofoil chord and hence on the Reynolds number. This can be 
overcome if the shear stresses on the tunnel floor and ceiling and the momentum 
thicknesses of the boundary layers entering the control surface are also measured, 
or if measurements are also taken along horizontal segments of traverse line in the 
flow, or along a curved traverse line altogether; in either case, equation (1) needs 
modification. During the development of the traverse method, however, the vertical 
traverse line crossing the nose shocks is likely to be adopted for the sake of 
simplicity, and it will be assumed in the following that the traverse line is normal to 
the undisturbed stream, and that the distance, D, between the leading edge and the 
traverse line is restricted, so that the traverse penetrates into uniform flow upstream 
of the nose shocks (Fig. 1). (It will be noted that for a precise estimate of the 
maximum admissible value of D, the upstream influence™ of the shocks on tunnel 
floor and ceiling must be taken into account). 


A second condition, suggesting a further restriction on the aerofoil chord in 
some cases, arises from a consideration of the wake structure. For blunt-based 
aerofoils, the mixing process near the base has an important influence both on the 
development of the wake and on the flow pattern outside the wake". When the 
aim of the traverse is a drag measurement, rather than an exploration of these flow 
processes near the trailing edge or base, it appears desirable to locate the traverse 
downstream of the region in which they take place. Moreover, the interpretation 
of wake traverse results in terms of drag is simplified if the classical boundary layer 
approximation applies to the wake at the traverse position (Section 5.1). It follows 
that the chord should be appreciably less than D for blunt-based aerofoils with 
initially laminar wake at relatively low test-section Mach numbers, but that little 
additional restriction on the chord is required for the more usual case of sharp-edged 
aerofoils at small incidence and with turbulent boundary layer. In either case the 
Reynolds number is restricted to roughly half the value admissible for experiments 
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for which only the pressure distribution on the aerofoil surface needs to be indepen: 
dent of the height of the test section. 


The remaining tunnel interference is due to two causes. First, the flow is not 
completely uniform even in the empty test section, and a correction for the back- 
ground flow pattern is usually desirable. Secondly, interference may arise from 
the interaction of the aerofoil shocks with the boundary layers on the side walls. The 
shocks are normal to these walls and no direct reflections of the shocks occur, but 
the local thickening of the boundary layers may be expected to induce “ ridges ” of 
pressure, each contained in the narrow space between two envelopes of Mach cones. 
The traverse line therefore passes upstream of all the pressure ridges, and no inter- 
ference occurs, if the width of the test section exceeds 2D/(M? - 1)'/? and, in any 
case, if the test section is wider than it is high. If not, false pressure readings will 
be caused only locally, where the traverse line crosses a pressure ridge, and it may 
be possible to choose D so that this does not occur where the accuracy of the read- 
ings is particularly important, or even to estimate the amount of interference from a 
comparison of traverse results obtained for two values of D and hence to correct the 
readings approximately. 
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3. Wave Drag Contribution 
3.1. DRAG FORMULA 


By continuity, f (p,U — pu) dy=0, and hence (1) may be written 


0+w 


= [p, -pt+pu (U —u)] dy, 
pu(U (U - (3) 


and, by (2), Cw=Cr= 2 (yM?)-! [1-2 


The gas will be assumed perfect, with constant specific heats, so that the 
equation of state, 


=(p/p,)” exp {= = ‘ ‘ (4 

Pp _ pT 
A (4a) 


holds throughout the flow. The local stagnation pressure p, is the pressure that 
would be obtained if an element of fluid was brought to rest from its local condition 
by an adiabatic and reversible process; from (4) and (4a) therefore. 


The energy equation T,=4q'? +c,T =constant ‘ (5) 

or, by (4), =constant, . . . . (5a) 


holds throughout the flow outside the wake. By (5) and (4) therefore, 


Pe { - , ~ « « 


outside the wake. Equations (4), (5a) and (5b) permit p and q to be expressed in 
terms of static and stagnation pressure, but not u in (3). 


Since aerofoils are hardly of practical importance unless their drag coefficients 
are small, it will now be assumed that the aerofoil is aerodynamically thin, with 
thickness parameter denoted by >. The perturbation, p’=(p—p,)/p,, of static 
pressure is then related to the perturbation 4 of stream direction by (Ref. 1, p. 121) 


with m?=M?-— 1, both across shocks and in the flow downstream of shocks outside 
the wake. 
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Expansion of (4), (Sa), (5b) and (6) in powers of the perturbations and substitu- 
tion in (3), noting that p’=O (r), but p,’=(p, po, =O (7°), yields 


It may be noted that the validity of (7) is not affected by the occurrence of an 
extensive interaction between tail shock and boundary layer, with marked bifurca- 
tion of the shock. In that case, each of the shocks near the aerofoil is even weaker 
than the joint shock farther out, and (6) therefore remains valid. 


On the other hand, if the nose of the aerofoil is blunt, the nose shock is detached 
and part of it is not weak. Unless the aerofoil has a practically excessive drag, 
however, the shock becomes weak at a certain distance from the nose, and there will 
be a pair of streamlines, one on each side of the aerofoil, and not far distant from it, 
beyond which (7) holds. The strips between these streamlines and the 
wake proper may be regarded as a special shock wave, in which | p,’ | is not <<< 1 
and which is to be included in the wake traverse. 


3.2. NUMERICAL COMPARISON 


It is useful to supplement (7) by a rough numerical estimate of the relative 
magnitude of the terms in that equation. The double-wedge aerofoil furnishes a 
convenient example, since for it the surface pressures are representative of the peak 
(static) pressures on the traverse*. A first approximation in terms of thickness ratio 
z and incidence z is obtained from (6); the first term on the right hand side of (7) is 


, 


and the second term, 
[24 (y— 1+ (y+ 1)/6) 
= —M? [2+(y—1+(+1)/6) (M? - 1)-°? (72), 
in the region of uniform flow downstream of the nose shock. Some values are given 
in Table I, for M =3 (the variation of w., with M is small for 2 <M < 4); the fourth 
and fifth lines of the table refer to the upper and lower side, respectively. It is seen 


that neglect of w,, corresponding to the neglect of 6? against unity in Refs. 1 (p. 562) 
and 2, implies a serious error. 


*It is assumed that, as in most cases, the traverse line crosses all the regions of uniform 
flow between the shocks. For aerofoils with curved surfaces, the peak pressures on the 
traverse depend on the traverse position, are lower than those on the surface and may often 
be expected to be better represented by those on the surface of the double-wedge aerofoil 
with the same thickness ratio and incidence. 
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TABLE I 
10°w,  —10°~, 10°w, 

for M=2 3 4 
7=0:06 
7 08 03 O1l O13 016 
7=0:10 
} 20 3704 018 O21 026 
7=0:10 0-04 
en 18 7 | 031 036 0-45 


The first approximation to the relation between entropy rise and shock deflec- 
tion (Ref. 1, p. 122), together with (5b), gives for the third term on the right hand 
side of (7), 


in the region of uniform flow downstream of the nose shock. Some values are again 
given in Table I (for M=3). It should be noted, however, that for curved aerofoils 
these values are indicative only of the lowest values of w, occurring on the traverse, 
since (- p,’) rises along it, from the nose shock to the tail shock. If the distance 
between trailing edge and traverse is small, the value of w, near the tail shock may be 
more nearly represented by the value for the double wedge, not of the same, thickness 
ratio, but of the same leading edge angle. Moreover, if the static pressure varies 
linearly with traverse distance, between nose and tail shocks, as it does approximately 
for curved aerofoils, the average of w, over this traverse segment is only a third of 
the figure given in Table I for <=0. A better, although still conservative, indication 
of the relative importance of w, over the traverse segment between the nose and 
tail shocks is therefore given by the minimum of the ratio 3w,/w,, of which some 
values are also shown in Table J. It is seen that the stagnation pressure distribution 
between nose and tail shocks makes a more important contribution to the drag than 
(7) might suggest and that the importance of this contribution rises markedly with 
thickness ratio and, especially, with incidence, and also to a lesser degree with Mach 
number. A closer approximation to cw than (7) therefore appears desirable in 
many cases*, 


Downstream of the tail shock, the third term in (7) becomes even more impor- 
tant than between the shocks. The static pressure perturbation, on the other hand, 
falls to very small values, and a satisfactory approximation to cw may then be 
expected to be given by equations (14) or (16) in Section 5.1 with T,’=0. 


4. Measurement of the Wave Drag Contribution 


4.1. STATIC PRESSURE 

A measurement of p - p, along the traverse is required for the determination of 
the first two terms in (7), and the major contribution to the integrals over the two 
terms is furnished by the largest values of | p-p,|. Hence the most serious effect 


*Note added in proof: Such an approximation may be obtained from the theory of Ref. 14. 


May 1957 131 


1) 
n 
Tt 
e 


R. E. MEYER 


arises from the measuring errors on the traverse segments where | p’ | is largest, 
i.e. from the segments near the shocks, and it is just on these segments that static 
pressure is not a measurable quantity, because of the interaction between shock and 
boundary layer on the probe! 


A conceivable approach would be to measure the pressure over most of the 
traverse, where it is measurable and known to be a smooth function of length along 
the traverse, and then to extrapolate this measured distribution to the shocks, having 
determined their position by a different, e.g. optical, method. It is instructive to 
begin the discussion by determining the order of accuracy that would then be 
required for the location of the shocks. 


For simplicity, assume that the pressure varies linearly with traverse length 
between the nose and tail shocks and that | p’| reaches approximately the same 
value at both shocks. If y” denotes traverse length measured outwards from the 
mid-point of this traverse segment, the length of which is approximately 
c (M? - 1)~'/?=c/m, the average of w, over the segment is 


c/Qm) 
w, | yay’. 
—c/(2m) 


A relative error ==2m A y”’/c in the ordinate of the nose shock (on both sides of 
the aerofdil) then causes a relative error 3:/2 in w,, approximately. A similar error 
may occur in the location of the tail shock, and since the two errors may add up, a 
total error 3< in w, must be anticipated. 


The importance of such an error depends, of course, on the relative contribution 
which w, makes to cp. On the other hand, sources of appreciable error will be seen 
to occur in most of the measurements required for the determination of the drag, 
and the most reasonable aim to bear in mind at the present stage appears to be a 
limitation of all of them to the same relative level. The relative error in cp will 
then also be 3: and, for instance, for M=3 and a chord of 3 in., an error 3| «| << 
5 per cent in cp implies an absolute error | 4 y” | < 10~? in. in shock location. It 
does not appear an easy matter to achieve such accuracy by optical means. 


Another method of shock location suggesting itself is the deduction of the shock 
position from an analysis of the distribution of static pressure recorded by the probe 
as it traverses across the shock. If a suitable probe is chosen, a flat-plate boundary 
layer and a two-dimensional interaction between it and the shock can be obtained. 
Whichever rule be chosen to deduce the approximate shock position from the 
recorded pressure distribution, the possible error will be some fraction of the total 
distance over which the interaction extends and, hence, some fraction € of the 
“* distance of upstream influence”, d, of the shock on the probe. The corresponding 
absolute error in the traverse direction is A y”’=éd/m and the relative error 
¢=2¢d/c approximately. A common case is that where the boundary layer suffers 
transition during the interaction, and then 


(Cyn), 
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where 5, is the displacement thickness just upstream of the interaction, R, the 
Reynolds number based on the probe length / upstream of interaction, c;, the 
pressure coefficient of the total pressure rise occurring during the interaction and 
F (cpm) a function plotted in Ref. 3. The displacement thickness is (Refs. 1, p. 416 
and 8, p. 136) 

), 


at least for e=1 and »/T=constant, and hence 
(1+0°277 M?) F (cpm). 


For a double wedge of thickness ratio ==0-1 at zero incidence, M=3 and R,.=10°, 
for example, with //c=1/3, say, ~0-28€ and thus 3|<| <5 per cent, as in the 
earlier example, only if it can indeed be ensured that £ < 1/16. 


The error in cp, however, may be expected to be smaller than 3< if the second 
method is used. If the strength of the tail shock is the same as that of the nose 
shock, then the error is of approximately the same magnitude, and in the same sense, 
for both shocks and, if it amounts to an overestimate of w, at the nose shock, then it 
amounts to a similar underestimate of w, at the tail shock and the total error 3 in w, 
is much less than 3<. For an aerofoil at incidence, the four shocks may be grouped 
into two pairs of approximately equal strength and a similar compensation occurs. 


It appears advantageous, therefore, to direct effort towards ensuring that the 
measuring error in static pressure near the shocks is not of a random, but of a 
systematic nature, so that compensation occurs, rather than towards reducing the 
individual errors to the desirable level. The method just proposed for achieving 
this may not be an efficient one, but it indicates that a useful level of accuracy should 
be attainable, even though compensation of the errors with regard to w, implies that 
they add up with regard to w,. 


Measuring errors in static pressure must, of course, be anticipated also where 
the static pressure is indeed a measurable quantity, but their effect on w, may be 
expected to be of less importance. On the other hand, it must not be forgotten that 
the undisturbed static pressure p, and the undisturbed Mach number M are also 
measured quantities. 


4.2. STAGNATION PRESSURE 


By (4b) and (5), the local stagnation pressure may be expressed in terms of local 
static pressure and local Mach number by 


The Pitot pressure which can actually be measured is the stagnation pressure down- 
stream of a normal shock at the local Mach number (Ref. 1, p. 556), 


Pr=Po {q+ (2yM? y+ (1+ 1) M?/2] 


= pf (M). 
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Consider first the section of traverse between the nose and tail shocks. Since 
=O and M,- M=O (7), =O (7) and the stagnation pressure perturbation 
is obtained from (9) as the small difference of two very much larger quantities; for a 
double wedge of s=0-1 at zero incidence and M=3, for example, p;’ ~ —30p,’ 
between the shocks. The Pitot tube appears to measure static pressure, between the 
shocks, rather than total head! 


The situation is similar to that encountered in (3), where cw (which is O (r*)) is 
represented in the form of a difference of quantities O (7). The second-order theory 
of the flow, which enabled (3) to be transformed into the more satisfactory form (7), 
does not, however, suffice to transform p,’ into a form better than a difference of 
quantities O (7), and further theoretical work is required. 


Nevertheless, the application of the traverse method need not await completion 
of such work in all cases. For a double wedge, p,’=constant between the nose and 
tail shocks and its value can thus be obtained directly from the measured pressure 
ratio of the nose shock. For curved aerofoils, p,’ can be measured at three 
points : —(i) just downstream of the nose shock, by the help of the measured static 
pressure ratio, (ii) where p’=0 between the shocks, by the help of a Pitot tube and 
(8), (9), and (iii) just upstream of the tail shock, by the help of the measured static 
pressure ratio and the measured stagnation pressure just downstream of the shock, 
provided that | p’| does fall to a very small value immediately downstream of the 
shock. Since p,’ falls smoothly along the traverse, from nose to tail shocks, a know- 
ledge of its values at these three points suffices to define its distribution, and 
although no great accuracy can be expected from the three measurements*, they may 
suffice for thin aerofoils at very small incidence, where the third term in (7) is only 
of minor importance and may thus be regarded as a correction term to be estimated, 
rather than measured. 


Similar difficulties will not, in general, be encountered on the section of traverse 
between tail shock and wake, since inviscid theory indicates very small values of p’ 
there, and a first approximation for cw in terms of py’ and p’ is derived in Section 5.2 
(equation (21))._ An exception occurs when the tail shock-boundary layer interaction 
is very extensive, so that p’ retains appreciable values even downstream of the tail 
shock and the approximation used in the following section becomes valid only near 
the edge of the wake. 


5. Friction Drag Contribution 


5.1. DRAG FORMULA 


Let u and wv denote the components of velocity parallel and normal to the 
undisturbed stream. The assumption will be made that v? may be neglected in 
comparison with u? throughout the wake. This assumption appears physically 
justified, provided that the traverse is not made too close to the trailing edge or base, 
for the flow outside the wake may be expected to approach a near-uniform flow 


*e.g. the error « in the measurement of p’ near the shocks discussed in the preceding section 
implies an error 3¢ in the value of p,’ there, since p, =O (p’). 
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rapidly with distance from the trailing edge, so that v?<< uv? at the edge of the 
wake. Some physical reasons which make it desirable to place the traverse at a 
sufficient distance from the trailing edge or base are indicated in Section 2. More- 
over, the wake will then be thin at the traverse position and the boundary layer 
approximation will apply, so that v? may be neglected against u? throughout the 
traverse of the wake. 


The stagnation temperature is therefore 


u? 


T,=T+ 2c,’ 


(10) 
and, by (4a), 


and, by (3), 


= 
4yM?cp= aT 


where B=(y-1)/y. 


It is convenient to note here that, also by (4a) and (10), 


_ (1+p’)(7,/T) (T,- (12) 
(Pu? 
By (4b), r=7,(2)" a+r.) 
(13) 
and T,-T=7,(1+7,[ 7" (2) 


Since the boundary layer approximation applies, the static pressure is constant 
across the wake, and since v? << u? at the edge, | p’ | << 1. The Prandtl numbers 
of gases do not differ much from unity, and hence | T,’ | << 1 in the wake, if there 
is little heat transfer across the aerofoil surface. If higher powers and products of 
p’ and T, are neglected, (11) to (13) yield 


cot =2 Dos)? (18) (1 + (7) +O +0 (T,”), 


(14) 
where g=1-—(1-A,/A,)'/?, 
BA, =(Po1/ - 1, 
BA, =T,,/T;, - 1=(y— 1) M?/2, . 
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By the energy-integral equation of boundary layer theory (Ref. 1, p. 398) 


J peau dy=O, 


w 


the heat transferred from the aerofoil to the fluid per second and unit span, and 
hence 


| puT , 


Ww 


(15) 


01 


where Ks is the Stanton number characterising the heat transfer from the aerofoil 
and 7,, is a representative wall temperature. In the absence of heat transfer, there- 
fore, the right hand side of (14) is equivalent to cp for the evaluation of the drag 
from (2). On the other hand, when heat transfer is not completely absent, but is 
small enough for | 7,’ | to remain small, neglect of stagnation temperature measure- 
ment will result in the first place, by (14) and (15), in measurement, not of the drag 
coefficient, but of the sum of drag coefficient and a certain fraction of the Stanton 
number. 


Equation (14) is also the formula for cp corresponding to the result of Ref. 9*, 
where tables will be found of A, against (p,, - p,)/Po=— Po’ /(1 + p,’) and of g against 
A,/Ao- 


Near the edge of the wake, | p,’ | << 1 and a more useful formula, 


, 2 , 
pul,” _ (M? 1) (p.'/2+P’) T, -| +0 (p,), (16) 


p,U yM? yM? 

is obtained by expanding (14). It is seen that the leading terms of (7) and (14) are 
the same near the edge of the wake. In fact, since it is based on the neglect of ~” 
against uw? only, and not on (6), the right hand side of (14), with T,’=0, may be 
expected to represent a better approximation than (7) for cw over the traverse 
segment from the edge of the wake to the neighbourhood of the tail shock. 


Figure 2 shows equations (14) and (16), for p’=7,’=0 and M=2-41, and it is 
seen that (16) provides a good approximation over an appreciable range of values 
of p,’. 

*The two results become identical when 7,’=p’=0, but the term O(p’) found here is not 
the same as that found in Ref. 9, where the coefficient » of p’ takes the value 
4{(v+1)/y-(1+8A,)/(,—A,)]._ It may be noted that, apart from the neglect of higher powers 
and products of p’ and T,’, the only assumption involved in the derivation of (14) is that 
v? may be neglected against u?. This assumption is also made in the subsonic case'®) where 
it also restricts the minimum distance between trailing edge and traverse, unless incidence 
and trailing edge angle are small. Further assumptions made in Ref. 9, but not here, are 
that all variations of stagnation temperature may be neglected and that the variation of 
entropy along mean streamlines in the wake may be neglected even over large distances. 
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5.2. PrroT EVALUATION 


It remains to express the stagnation pressure in terms of the Pitot pressure, 
defined by (9), and the static pressure. By (8), which is now based on (10), rather 


(17a) 


than (5), 
1+(y- 1) M?/2=(14+(- 1) M?/2)(1+p,. (1+ p)-8 (17) 
so that M.=M, (p,’, p’; M) 
OM, _ 
and (1+ Op.’ = 


It is convenient to introduce an auxiliary Mach number, 
M,=M, (p,’,0; M) 
=M (1-8) 
by (17) and (14a). Then, by (17a) and (14a), 


dM 
M.- - p’ (1 +p.) ap,’ 


[1+ (y- 1) M,?/2] 
yM, 


to the first order, if | p’ | <<1, 
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Pr _ Po 
Finally, by (9) and (18a), Pr, Po. f(M) 


=h(p,.’;M)[1- p’ j 


_ Po {(Mn) 
where /A(p,’; M)= rw 7(M) 


and 


= —2(M,? - 1)?M,-? (2yM,? - y +1)" 


by (9). Fig. 3 shows h and 1+/j against p,’ for M=2-41. 


(19) 


(19a) 


The best procedure for obtaining p,’ from the recorded values of pr’ and p’ 
would appear to be as follows. First, the functions M, and h are tabulated against 
p,’, for the Mach number M of the nozzle, by the help of (18), (14a), the tables of 
Ref. 9, (19a) and a table of f (M) (given, e.g., in Ref. 10 as Pyn/P,.). Secondly, an 


auxiliary quantity, p,”, is determined from 


Pr, =h (p."; M), 
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and finally, Po’ + kp’ +O (p””), 


” ” d 
where k=] (py ;M)h(p, ; M) 
0 


= (1+ p,”) (M,?— 1)? (1+ 


with M,, understood as M, (p,”; M), by (19), (19a) and since 


dh(p,’) _ {(M,) f (Ms) 


by (19a) and (18a). The values of p,/p,, for which M,=1 are 0:24, 0-052 and 0-012 
for M=2, 3 and 4, respectively, and subsonic velocities will therefore hardly be 
encountered at the traverse position. 


By (20), h=1+(1+j.) po.’ +O (p,”), 
near the edge of the wake, with 
jo=j (0; M) = —2 (M? - 1)°M-? (2yM? - y+ by (19a), 


and the first approximation for cp and cw near the edge of the wake is by (19a), 


T 0 T e 21 
+ jop') (1 +0 + O +O (T.), (21) 


6. Measurement of the Friction Drag Contribution 


6.1. STAGNATION TEMPERATURE 


For fluids with Prandtl number o ~ 1, little seems known to date about the 
temperature distribution in a wake due to dissipation, in the absence of heat transfer. 
It appears plausible, however, that the main process in the wake is a decay of the 
velocity and temperature gradients produced by the aerofoil, and the excess and 
deficiency of stagnation temperature found near the aerofoil may therefore be 
expected to provide bounds for the excess and deficiency in the wake. The results 
for the flat plate, in turn, may plausibly be expected to provide a guide to the order 
of magnitude of 7,’ for thin aerofoils. Experiment“” indicates, moreover, that the 
recovery factor, (T~—7,)/(T,,—T,), is higher for the turbulent boundary layer on 
the flat plate than for the laminar layer, and the largest values of | 7,’ | may there- 
fore be expected to be those encountered with the laminar layer, where (Ref. 1, 
p. 418) 

1) M? 
2+(y-1)M? 
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the value occurring at the wall. For «=0-72, therefore, T,’”/8 <0:003 near the 
flat plate and measurement of this term in (14) appears unwarranted, in the absence 
of heat transfer. 


Since special arrangements to avoid heat transfer are not usually made, its 
complete absence is however not always assured. The most common case for 
experiments in intermittent tunnels is that the stagnation temperature and the initial 
temperature of the model are both atmospheric. Since the “Thermometer” 
temperature at the model surface due to dissipation is then lower, a certain amount 
of heat transfer from the model to the air must take place*, at least during an early 
part of the run. To estimate how rapidly the rate of heat transfer falls during the 
run, assume that the temperature distribution in the model is uniform at any instant 
but, for simplicity, neglect the heat conducted into it from the support. The heat 
transferred to the fluid per second and per unit span is then 


aT. 
Q= —CmP nV 

where Cm, Pm and 7, are respectively the specific heat per unit mass, density and 
temperature of the model and V=Arc? is its volume per unit span. The Stanton 
number, Ks, of the aerofoil is defined by 


Q 
cp p,Uc(T.-T,)’ 


and since it is proportional (Ref. 1, p. 426) to (T~— Tn)/(Tw—T,), where T,, denotes 
the thermometer temperature, 


Tx Tin 
(22) 


p,UcKs’ 


where k= ota’ 


CmPm Por Pa ATC "7 


(T~-T,) 


(23) 


*Such heat transfer implies, of course, that the deficiency of stagnation temperature will be 
even smaller than it is in the absence of heat transfer. 


The stagnation temperature T,, is assumed constant. For blow-down tunnels without tempera- 
ture control, the right hand side must be multiplied by 
t 
1—(2+r (y—1) M?) (2+(y—-1) e* / dt) dt, 
0 
where r is the recovery factor and T,,, (¢) the stagnation temperature in the test section. 
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p, and p, denote respectively the atmospheric density and pressure of the air, U the 
absolute velocity in the test section and c,p,/(Cmpm)=3°5 x 10~* for air and steel. 
The results for the flat plate may be expected to provide an estimate for the Stanton 
number also for thin aerofoils. For the double-sided flat plate with laminar 
boundary layer (Ref. 1, p. 426), 


where C; is the average skin-friction coefficient for one side and R, the Reynolds 
number based on the chord. For the turbulent compressible boundary layer Ks’ 
cannot yet be predicted with similar certainty, but if the low-speed results provide 
a guide in a similar way as they do in the laminar case, a sufficient approximation 
for the present purpose is obtained by assuming 


As an example, consider a double wedge (A=1/2) of 3 in. chord and thickness 
ratio r=0-1 in a stream of Mach number 2:4. In a suck-down tunnel (p,,=p, and 
R.=0-75 x 10°, approximately) k=0-06 for a completely laminar boundary layer, 
and k=0-1 for a completely turbulent layer; in a blow-down tunnel with p,.=4p, 
(and R.=3 x 10°, approximately) the corresponding values are 0:03 and 0-075. In 
either case, therefore, the temperature difference T,,— 7 falls to half its initial value 
in a time of 10 to 15 seconds and the process of heat transfer from model to stream 
is neither virtually completed in the first few seconds of the run, nor sensibly 
constant during the run. 


It follows firstly that the drag differs from that at zero heat transfer, which is 
intended to be measured, and secondly that an error is introduced into the traverse 
method by neglect to determine the second term on the left hand side of (14). An 
estimate of the first error can be obtained, for the case of a laminar layer with 
»/T°=constant, from the approximate formula for the flat plate (Ref. 1, p. 422), 


C,=1-328R.-*/? [0-45 + 0:55T /T, +0-09 (y 1) 


Since (Ref. 1, p. 418), 
Ta/T,= 1+(y-1)0°/?M?/2 


and by (14a) and (22), the relative error in C; due to the heat transfer is 


0:275 (1-0?) (y— 1) 
2+0°730%!? (y—1) M? 


For the turbulent layer, a similarly confident estimate does not appear possible at 
present, but a relation that has been suggested, on the basis of experimental 
evidence, 

C:=(T,/Tw) Cu “Rel, 
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where C;[R;] denotes the incompressible, average skin-friction coefficient for the 
flat plate of Reynolds number R;, which may be approximated (Ref. 8, p. 362) by 
C;R;'/*=constant. If the thermometer temperature is taken to be (Ref. 1, p. 456) 


this leads to a relative error 


= M? 


in C, due to the heat transfer. For example, for »w=8/9, c=0°72 and M=2-4, 
€,= —0:003le-*' and ¢,= —0-025e-**. 


The second error, which consists in measuring Cp + Ks’ (Ty~—T1)/T,, instead of 
Cr, by (14), (15) and (23), can be similarly estimated. For the flat plate with laminar 
boundary layer, the relative error in Cy is 


M) 


(28) 


by (24), (26), (22) and (14a). For the turbulent layer, if (25) and (27) are taken 
instead of (24) and (26), the error is 


(1-0!) (y— 1) 
2+(y-1)M? 


Thus for «=0:72 and M=2:4, for comparison with the preceding example, 
8,=0-050 e~*' and 6,=0-028 e-**. It is seen that the two errors introduced by the 
heat transfer always counteract each other and, for mainly turbulent boundary 
layers, may indeed be expected to compensate each other effectively. For pre- 
dominantly laminar layers, the second error is a significant one, and an approximate 
correction of the results by the help of (28) and (22) appears desirable. 


6.2. PiroT PRESSURE 


The curve of cy against y is roughly of Gaussian type, i.e. 
cry ~aexp(- b*y?/c?) 


and dy/c =~ at'!/?/b, 


w 


and errors in the recorded pressures may lead to errors in a, i.e. in the height of the 
curve, or in b, i.e. in the half-width of the curve, or in both. 


Assume first that the half-width is correctly indicated by the Pitot readings. A 
relative error 6 in a, and hence also in Cy, then corresponds to an error ¢«=5/¥ in the 
recorded values of py’, where 
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— Pr 

Cp Opy 

hel (2g—1) (1+ 


by (19), (20), (14) and (14a), if T,’ and p’ are neglected. The magnification factor v 
remains close to unity in the practical range of stagnation pressures (Fig. 3) and 
thus, if the readings indicate a Pitot pressure distribution in excess of the correct one 
by a fairly uniform percentage, the total friction drag contribution will be over- 
estimated in approximately the same proportion. 


Errors in a may arise, e.g. from small day-to-day changes in tunnel Mach 
number, and differences of only +0-005 in M were found in Ref. 3 to lead to differ- 
ences of up to 5 per cent in a. Close repetition is therefore important for traverse 
measurements. 


On the other hand, an error of the “‘ displacement” type®.’® due to the trans- 
verse gradient of total head will affect primarily the half-width of the friction drag 
contribution. It should be noted that the measurement of the friction drag by Pitot 
traverse is very sensitive to this type of error. Since v~1, the half-widths of 
friction drag contribution: and Pitot pressure distribution are approximately the 
same and hence, if the readings indicate half the maximum value of (— py’) at a 
traverse distance from the centre of the wake which is in excess of the correct value, 
(c/b) (log 2)'/?, by a fraction 5, then Cp will be over-estimated in the same proportion. 
Since the wake is narrow, the absolute displacement permissible for adequate 
accuracy of Cy is therefore extremely small. It would appear that further experi- 
mental work on the behaviour of Pitot tubes in narrow supersonic wakes is desirable 
before adequate accuracy can be expected confidently from a traverse measurement 
of the friction drag of thin aerofoils at small incidence. 


A relative error ¢ in p’, finally, leads to a relative error 5 in cp such that 


p’ 


approximately, if | p’|<<1, by (14), (14a), (19), (19a), (20) and (29). Thus 
5=O (ep’/ py’) near the edge of the wake, but only O (<p’) in the main part of the 
wake, and Cy is not sensitive to errors in p’. Indeed, if | p’ | is as small at the edge 
of the wake as inviscid theory suggests, a static pressure traverse of the wake region 
appears unnecessary. 
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A Reduction Method for Buckling Problems 
of Orthotropic Plates 


P. SHULESHKO*, A.M.I.E.Aust. 


(New South Wales University of Technology) 


SumMMaRY: Several plate buckling problems are solved, using a reduction 

method. By this method the solution of an orthotropic plate can be reduced 

to the solution of an isotropic plate and the solution of a plate with bi-axial 

loading can be reduced to the solution of a plate with uni-axial loading and 

so on. Plates with simply-supported ends and various boundary conditions 

at the sides with uni-axial and bi-axial loading are considered and the necessary 
reduction equations are given. 


1. Introduction 


The present paper presents an application of the reduction method to the 
solution of problems of the buckling of rectangular orthotropic plates under 
combined tension and compression. 


This method was developed by the author in 1954 at the N.S.W. University of 
Technology and a paper on it was read to the Shevchenko Scientific Society 
(Proceedings No. 4, 7th December, 1954). The method is based on the use of 
simple simultaneous equations which are called reduction equations. By using 
these equations a complex problem can be reduced to a simple problem. In this 
paper the problem of finding the critical load for an orthotropic plate is reduced to a 
problem of finding the critical load for an isotropic plate. 


NOTATION 
a,b _ length and width of rectangular isotropic plate 
a,,b, length and width of rectangular orthotropic plate 


x,,y, rectangular co-ordinates 


See Fig. 1 
x,y = x,/a, and y,/b, respectively 


D _ bending stiffness of the isotropic plate 
D,,,D,,,D,, bending stiffnesses of the orthotropic plate 


*Formerly Professor of Theoretical and Applied Mechanics, Ukrainian Technical University, 
Germany. 
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w lateral deflection of a point in middle surface of the buckled 
plate 


P,,P, compressive forces per unit length applied in the x- and 
y-directions 


Y _ function of y alone in equations (5) and (6) 
C,,C,,C;,,C, constants of integration in equation (7) 
#  Poisson’s ratio of the contraction of the isotropic plate 


#,  Poisson’s ratio of the contraction of the orthotropic plate in the 
y-direction to the extension in the x-direction resulting from pure 
tension in the x-direction 


S, stiffness per unit length of the elastic restraining medium at the 
edge of the isotropic plate, or moment required to rotate a unit 
length of the medium through one-fourth of radian 


S, stiffness per unit length of the elastic restraining medium at the 
edge of the orthotropic plate, or moment required to rotate a 
unit length of the medium through one-fourth of radian 


m,m, number of half waves of the buckled isotropic and orthotropic 
plate, respectively, in the x-direction. 


Orthotropic Plate Isotropic Plate 
a, B (see relations (8)) 2, 8, (see relations (12) ) 
a,=D,,//(D,,D.2) a= 
a,=(D,,/D,,) a, = 
A, = me = me 
P.a,’ _P,a?_ 
WFD 
r,=4S,b,/D,, r,=4S,b/D. 


2. The Differential Equation for an Orthotropic Plate 


Assume that forces P, are distributed along the edges x,=0 and x,=a, of the 
plate and that forces P, are distributed along the sides y, =0 and y, =), (Fig. 1 (@) ). 
Forces P, and P, can produce a compression or tension of the plate. 
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Fig. 1. 


Rectangular plate compressed in two perpendicular directions. (a) Original plate. 
(b) Transformed plate. 


A differential equation of the orthotropic plate can be written, as a particular 
case of a more general equation given in Refs. 1, 2, 3 and 4 in the following form. 


4 : 4 2 2 


+2D, +D,.—— +P 


where D,,, D,. and D,, are the bending stiffnesses of the plate." * > * ” 


Equation (1) can be transformed into non-dimensional form by using the 
substitution 


After transformation, the edges of the plate x,=—0, x,=a, and sides y,=0, 
y,=b, will be changed to x=0, x=1 and to y=0, y=1 (Fig. 1 (6)). The 
differential equation (1) becomes 


1 


Di; = (Pu), Mb 
P,b,? 


(4b) 


k, k, = 


The problem is to find the least value of k,, for a given k,, (or to find k,, for 
a given value of k,, or for a given ratio v,=k,,/k,.) and A, for which equation (3) 
has a solution when w is not zero and which satisfies the boundary conditions for 


the given problem. 
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The solution will be different for different boundary conditions. Consider a 
plate with edges x=0 and x=1 simply-supported, but for which the sides y=0 and 
y=1 can have a different type of support. 


The boundary conditions at the edges x=0 and x=1 will be satisfied if the 
solution is represented in the form‘: * ?” 


w=Y (y) sin ‘ ‘ ‘ (5) 
where 


Y(y) isa function of y alone, to be determined later 
m, _ is the number of half-waves of the buckled plate in the x-direction. 


Substituting (5) in (3), the following differential equation for determining the 
unknown function Y is obtained: — 


1 


The solution of equation (6) can be written 


Y =C, sinh zy+C, cosh zy+C, sin By+C, cos By ‘ ; (7) 


C,, C,, C,, and C, are constants of integration which can be determined from the 
boundary conditions for the sides y=0 and y=1. 


3. Particular Cases of the Buckling of Plates 


Several cases are considered with different boundary conditions for the sides 
y=0 and y=1. 


3.1. THE SipeSs y=0 AND y=1 ARE CLAMPED 


The boundary conditions for sides y=0 and y=1 are 
and Yy’=0 ‘ ‘ 


After satisfying the boundary conditions (9) the following transcendental 
equation is obtained : — 


(cos 8B—cosh a)? + (sin B- sinh x) ( sin B+ =0 . (10a) 


or . 4 . (10d) 
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The transcendental equation for an isotropic plate“ will have the same form 
as equations (10) 


f 8.) =90, ‘ (il) 


but the expressions for z, and 8, will be 


(12) 
It is obvious that, if 


equation (10b) will be transformed into equation (11) and first problem will be 
reduced to the second or vice versa. 


Equations (13) will be fulfilled if the following simultaneous equations are 
satisfied : — 


Ay? 24,)=Ay? (Kos —2) (14) 


A,‘ (k,,A,-?- 1)=A,* 1) 
Equation (14) represents the relationships between roots of the two trans- 


cendental equations (10a) and (11). If the roots of equation (11) are known the 
roots of equation (10a) can be calculated by the use of equations (14) and vice versa. 


Taking A,=A,, then, from (14), 


+2 (@,—1) 


and vice versa. 
3.2. SOME PARTICULAR CASES OF EQuaATIONS (14) 
3.2.1. Uni-axial loading in y-direction only 


If the isotropic plate has uni-axial loading in the y-direction only, i.e. k,,=0, 
then, from (14), 


Ay=A, (1—k, (16) 


+ (ko.—2) kA. 
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Equations similar to equations (16) were obtained by another method by 
W. H. Wittrick®. 


As a rule, the values of A, and k,, in (16) must be given. The value of A, can 
be calculated from the first equation of (16). Then the value of k,,, for a calculated 
value of A,, can be obtained from a table or graph prepared for the isotropic plate; 
finally, the value of k,, can be calculated from the second equation of (16). 


Equations (16) can be used when k,, < A,? including all negative values of k,, 
(i.e. tensile). 
3.2.2. Uni-axial loading in x-direction only 


If the isotropic plate has uni-axial loading in the x-direction only, i.e. k,.=0, 
then, frorn (14), 


Ay=A, (a, --0°5k,.) 


1 
k,,=A,? +(a,—0°5k,.) (ko, — 


This case is possible when k,,<22,, including all negative values of k,. 
(i.e. tensile). 


Example 

Find a critical force for the orthotropic plate compressed in the x-direction if 
(b,/a,) =2, z,=0-40 and k,,=0. 
Solution 

The value of A, for the isotropic plate, from the first equation of (17) taking 
m,=1, is 
A, =2 x 1/(0-40-—0)= 1-265, 


1 


= 1965 =0-79. 


or 


From a table or graph, prepared for the isotropic plate compressed in the 
x-direction only (Fig. 2) 


k,, = 727. 
Then, from the second equation of (17), 
k,,= (2)? + (0-40—0) [7-27 —(1:265)?]=6-27. . . (18) 


The critical force, from the first of equations (45) is 


Py 6-27 On? 22) 
1 
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8 
7 
0-2 O-4 
| 
0-4 ole 0:8 
0-6 0-2 12 
+ 
0-8 1-2 1-6 


Fig. 2. 
Curve for the determination of the parameter k,, of the isotropic plate with simply-supported 
ends and clamped sides compressed in the horizontal direction by uniformly distributed forces 
P, per unit length of the end. 


Taking m,=2, then A,=4, A,=4/(0-40—0)=- 2°53. k,, =9°60 


and k,,=(4)? +(0-40— 0) [9-60 — (2-53)*]— 17-28. 


Then — 


(20) 


Hence the smallest value of P,.- is obtained by taking m,=1, i.e. by assuming 
that the buckled plate has only one half-wave deflection in the x-direction. 


The same value of k,,=6:27 was obtained by R. C. T. Smith (see his 
Table 1) from the solution of the complicated eigenvalue equation (11) of 
his paper. 

All the values of k,, (his A/=*) in his Table I for different values of A, (his 
B=m,/A,) and 2, (his z) can be obtained by using our equations (17). Part of 
Table I is given here, and the values of k,, are given side by side with Smith’s 
values of A/=* for comparison. 
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TABLE I 
1:0 | 11:92 11:92 | 7:69 7:69 8-61 8-61 
0:8 11°46 11°46 7:22 7-22 | 8:12 
06 | 11-00 10-99 6°75 6748 7:62 7-62 
0-4 10°54 1052, 627 6:27 7:13 
0-2 10-06 1006 | 5-78 5-78 6°64 


Table I shows very good agreement between the two sets of results. Using 
equations (17) it is not necessary to solve for every z, and 2, of the complicated 
transcendental equation (10a) or eigenvalue equation (11) in Smith’s paper. It is 
enough to solve equation (10a) or (11) only for «,=2z,=1 (isotropic plate) and 
k,.=0 (uni-axial loading). 


The derived reduction equations (16) and (17) and the example given show that, 
by using the Reduction Method, it is possible to get a solution for an orthotropic 
plate with a certain bi-axial loading if the solution for an isotropic plate with 
uni-axial loading is known. The same can be said about two isotropic plates with 
bi-axial and uni-axial loading. These important facts lead to a great reduction in 
the labour of computation. 


3.3. THE SIDE y=0 IS CLAMPED, THE SIDE y=1 IS SIMPLY-SUPPORTED 


The transcendental equation for this case''*: *) has the form 


and the reduction equations are 
a=a, and 
These equations will be fulfilled if the previous equations (14) are satisfied. 


Hence, for the three cases considered, the same analysis and the same equations 
(14), (16), and (17) can be used, except that the values of k,, and k,, have to be 
taken for corresponding isotropic plates with corresponding boundary conditions 
for the sides y=0 and y=1. 


3.4. THE Sipe y=0 IS FREE, THE SIDE y= 1 IS SIMPLY-SUPPORTED 


After satisfying the boundary conditions at y=0 and v=1 the following 
transcendental equation is obtained“ : — 


where B= 8? +A,7n, (23) 
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#, iS Poisson’s ratio of the contraction in the y-direction to the extension in the 
x-direction resulting from pure tension in the x-direction, and z and 8 are as given 
in equations (8). 

The expressions for A and B for the isotropic plate are 


A, = —A, ur’, B, B,? (24) 


The transcendental equation (22) can be transformed into the transcendental 
equation for the isotropic plate if the following equations are satisfied 


a=%,, . . 5) 
and A=A,, . . (26) 


Equation (25) will be fulfilled if the previous equations (14) are satisfied. 
Equations (26) will be fulfilled if the following equation is satisfied : — 


From the equation (27) and equations (14) the following expressions can be 


obtained : — 
a, 
ker = (ky, —A,”) (4 & 
The values of z,, z., “,. “, A, and k,, must be given for the actual plate; the 


values of A, and k,, can be obtained from equations (28) and (29). Then the value 
of k,, can be taken (for calculated A,) from the table or graph prepared for the 
isotropic plate®’. Finally, the value of k,. can be calculated from equation (30). 


3.5. THE SIDE y=0 IS FREE, THE SIDE y=! IS CLAMPED 
The transcendental equation for this case is 


sinh z sin B— 28 cosh z cos 8—228=—0 (31) 


where z, 8, A and B are as given in relations (8) and (23). The reduction equations 
are the same as those in Section 3.4. 
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3.6 THE SIDES y=0 AND y=1 ARE FREE AND UNIFORMLY LoaDED By Forces P, 


The transcendental equation for the case will have the form“ 


(32) 


BA®* sin 8 cosh z— zB" sinh z cos B=0 
or 8A? sinh z cos sin cosh a=0 


where 2, 8, A, and B are as given in relations (8) and (23). 


The first equation of (32) corresponds to the symmetrical forms of the buckled 
plate, the second to the anti-symmetrical forms. The reduction equations for this 
case are the same as in Section 3.4 (equations (28), (29), and (30) ). 


3.7. THE Sipe ¥y=0 IS FREE, THE SIDE y=1 IS ELASTICALLY RESTRAINED 
AGAINST ROTATION 


After satisfying the boundary conditions at the sides v=0 and y=1 the 
following transcendental equation is obtained“ : — 


(z | Sin 8 cosh x Be sinh sin | -r,| B aB cosh z cos 


- ( a? 4 — *) sinh z sin 8—226=0 (33) 


S, is the stiffness per unit length of the elastic restraining medium at the 
edge of the plate, or the moment required to rotate a unit length of the medium 
through one-fourth of radian (see Ref. 5). 


The reduction equations for this case are 


A=A,, B=B, ‘ . (36) 


These equations will be fulfilled if the previous equations (28), (29), (30) and 
the additional equation 


are satisfied. 
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3.8. THE y=0 IS SUPPORTED, THE SIDE y= 1 IS ELASTICALLY RESTRAINED 
AGAINST ROTATION 


The transcendental equation for this case is 
— + sinh z sin 8 +r, (6 sinh cos 8 —a cosh z sin 8)— 0 (39) 
The reduction equations are (14) and (38). 


3.9. THE SipE y=0 IS CLAMPED, THE SIDE y=1 IS ELASTICALLY RESTRAINED 
AGAINST ROTATION 
The transcendental equation is 


2? + B? 
ap 


(8 sinh z cos 2 cosh 2 sin 


[ (cos 8B —cosh +( sin sinh x) (sin B+ gsinh =0. (40) 


Equations (10) and (21) can be obtained from equation (40), as particular cases for 
r,=0O and r,=0. 


The reduction equations are the same as in Section 3.8. 


In the same way the reduction equations can be obtained for other cases. 


4. Conclusions 


From the examples of the buckling of plates considered it follows that, if the 
solution for an isotropic plate with simply-supported ends and with any type of 
support of the sides is known, the solution for an orthotropic plate with the same 
edge support can be obtained without difficulty, using the reduction equations 
given here. 


Sometimes it is possible to obtain a solution for a plate with a certain bi-axial 
loading if the solution for a plate with uni-axial loading is known (Sections 3.1, 
3.2 and 3.3). 


Cases where the ends of the plate have another type of support and the sides 
have any type of support will be considered in a later paper. 


The method considered leads to a great reduction in the labour of computation 
because it is not necessary to solve complicated transcendental equations or infinite 
determinants. The Method of Reduction can be applied not only to the theory of 
the buckling of plates but to many other problems in theoretical and applied 
mechanics and physics. 
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Matrix Inversion by Partitioning 
ERYK KOSKO 
(Avro Aircraft Ltd., Toronto) 


Summary: A systematic discussion of partitioning as a tool for matrix inversion 
is presented, together with various methods and applications which have been 
of help in actual computations. New concepts are introduced, among them 
those of super-matrix and of square partitioning. The most usual type of 
partitioning, that into 2 x2 sub-matrices, is discussed in detail, showing the 
orderly arrangement of the calculations in an auxiliary matrix. Further 
sections deal with matrices of the continuant type, and with special types of 
symmetry in the arrangement of the sub-matrices. The greatest advantage of 
the method of partitioning for the inversion of these types (as compared with 
the elimination method) lies in a considerable reduction in the number of 
arithmetical operations. 


1. Introduction 


The numerical inversion of large matrices is a problem which has important 
applications in many fields of science and technology. It usually involves laborious 
and costly calculations and much effort has been devoted in recent years to the 
development of time- and labour-saving methods. Among the methods available 
for inverting matrices, those based on partitioning seem to have received relatively 
little attention. A survey of texts on matrices"’"® reveals that partitioning 
and basic rules of addition and multiplication of partitioned matrices are mentioned 
in most of them, but only a few of the authors recommend partitioning as a method 
of inversion. Among those who discuss partitioning, Ref. 2 is the only one 
which demonstrates in a practical way the simplest rules. A search of periodic 
publications and reports shows that there is no systematic treatment of the subject 
in the literature. 


There appear to be three main cases in which one would expect partitionjng to 
be of some advantage when inverting a matrix : — 


(a) When certain blocks of elements are repeated within the matrix in a 
regular pattern or arranged in a symmetrical way, partitioning into sub-matrices 
comes quite naturally, as it merely follows the natural structure of the matrix (which 
is probably a reflection of the symmetric or periodic structure of the physical 
system represented by the matrix). The inverse matrix may also be expected to 
possess similar structural properties. 
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(b) Certain blocks of elements within the matrix may be of a special nature, 
e.g. they may consist entirely of zeros, or have the form of a unit sub-matrix or 
diagonal sub-matrix, etc. Operations with such sub-matrices may simplify the 
calculations. 


(c) Resort is often taken to partitioning when the order of the given matrix 
exceeds the capacity of the computing equipment. 


Much too often partitioning is used only for this last reason, without realising 
that the matrix considered belongs to some special category and that an appreciable 
number of arithmetical operations could be saved. The widespread use of large- 
scale computers appears to have developed a “ brute force” approach to problems 
which often can be tackled with desk calculators at much less expense. 


Similar lines of thought may apply to the calculation of characteristic roots and 
vectors, but this paper will be limited to the problem of inversion only. 


An attempt is made to present a systematic discussion of partitioning as a tool 
of matrix inversion. The most popular type of partitioning, in 2 x 2 sub-matrices, 
is discussed at some length. The calculations are best arranged in an orderly 
scheme called the auxiliary matrix. A similar scheme, based on recurrence 
relations, is demonstrated for the practically important matrices of the continuant 
type. Other kinds of symmetry are also treated, with a digression on the number 
of arithmetical operations required for inverting a matrix. 


As shown in Section 12, all the methods of inversion by partitioning referred 
to in the paper have one thing in common. They are based on a transformation of 
the given matrix into a diagonal super-matrix, the inversion of which amounts to 
the inversion of its diagonal blocks. The inverse of the given matrix is obtained 
by a transposed transformation of the inverse of the diagonal super-matrix. The 
inversion of a large scale matrix is thus reduced to the inversion of several smaller 
matrices (the sum of their orders being equal to the order of the original matrix). 
plus a number of matrix multiplications and additions representative of the trans- 
formations. Considerations of convergence, stability and accuracy affect mainly 
the matrix inversions and are therefore reduced in scale. 


All the derivations are elementary although quite general; they are based on 
the theorem of the uniqueness of the inverse of a non-singular matrix. The formulae 
presented can easily be verified by direct substitution and no further proofs 
are required. 


2. Definitions and Nomenclature 


For convenience, a square matrix partitioned into a number of sub-matrices is 
called a super-matrix, to emphasise that the sub-matrices, rather than the scalar 
elements, are to be considered as its constituents. The super-matrix may then be 
regarded as consisting of a number of super-rows or super-columns, a terminology 
which needs no further explanation. Some of the constituent sub-matrices may be 
single row or column vectors, or even scalar terms; some may be composed entirely 
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of zero elements in which case they are denoted by the usual symbol O. Scalar 
quantities are denoted by lower case letters, sub-matrices and ordinary (i.e. scalar) 
matrices by capitals, and super-matrices by bold face capitals. 


In this paper, the partitioning is always such that the sub-matrices containing 
elements of the principal diagonal (they are called here diagonal blocks) are square, 
so that their principal diagonal coincides with a segment of the principal diagonal 
of the super-matrix. Such partitioning will be called square partitioning. We may 
say then that the diagonal blocks form the super-diagonal of the matrix. If m is 
the number of diagonal blocks (also the number of subdivisions of any row or 
column) of the given matrix, then the total number of sub-matrices is m?, and we 
shall refer to ‘‘m x m’” partitioning or a partitioning of the m™ order. If moreover all 
the sub-matrices are square (and necessarily of the same order n, such that the order 
of the super-matrix is N=n x m), then we shall refer to a uniform partitioning. 


A simple application of these concepts is afforded by the diagonal super-matrix, 


a-[2 2], 


to 


sometimes also called the direct sum of the matrices B and C. In the ordinary 
(scalar) sense this cannot be considered, in general, as a diagonal matrix, for it may 
have non-zero elements outside the principal diagonal. But its inverse can be 
written as the direct sum of the inverse of its diagonal blocks 


in analogy to a well known property of diagonal matrices. Direct sums of more 
than two square matrices may be inverted in the same manner. 


Another example is that of a triangular super-matrix whose sub-matrices on 
one side of the super-diagonal are all zero. This is in general different from a 
scalar triangular matrix, since the diagonal blocks may contain non-zero elements 
on the null side of the principal diagonal. Such a triangular super-matrix can be 
inverted step by step, by successive inversions and substitutions, much like an 
ordinary triangular matrix. Clearly, the inverse of a lower (upper) triangular 
super-matrix is another lower (upper) triangular super-matrix. 


Let N be the order of the given matrix, and n,, m., ..., Mm, the orders of the 
m successive sub-matrices in the super-diagonal (or in any super-row or super- 
column). Then the m positive integers n; (i=1, 2, ..., m) must satisfy the relation 


+Mm=N. 


For the given matrix we may write A=[A,], (i, k=1, 2, ..., m), where the 
typical sub-matrix Aj, is of the order n; by n,. This representation applies only to 
Square partitioning, for the diagonal blocks A,; are necessarily square matrices of 
order n;. For uniformly square partitioning, all n,;=n=N/m. 
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As usual, the first subscript indicates the super-row, and the second the super- 
column, in which the sub-matrix is situated. 


With these conventions, there will be no need to refer explicitly to the order 
of the various matrices and sub-matrices. 


It is further assumed that the matrices which are to be inverted are all non- 
singular, and that this applies also to sub-matrices whose inverses may appear in 
the calculations. 


If a second matrix B of the same order as A is also partitioned m x m so that 
each sub-matrix B,, is of the same order n; x mn, as the sub-matrix A;,, the partition- 
ing of B will be referred to as corresponding to that of A. It is assumed that the 
inverse matrix P=A~' is always partitioned in a way corresponding to A. 


3. General Relations 


In matrix notation, the fact that P is the inverse of A is expressed by either of 
the equations 


AP=I or PA=I ‘ (1) 


where I is the unit matrix of order N. The first of these equations is equivalent to 
m? relations between sub-matrices 


l= 


each obtained by the multiplication of the i** super-row of A by the k" super- 
column of P. The alternative equation (1) is equivalent to another set of 
m? relations 


obtained similarly by multiplying the i" super-row of P by the k super-column of 
A. The symbol A; is a generalisation of the Kronecker 6 and represents either a 
unit matrix of order n; when i=k, or a zero matrix of order n;xm, when ik, 
The typical equation of the set (2) will be symbolised for brevity by (ix k’), and 
similarly an equation of the set (3) by (i’ xk), the prime (’) referring here to the 
inverse matrix. 


It is seen that any of the partial products AyP;, or PyAy, is a matrix of order 
(1; x m) x (nm, x ny) =(n; x n,); hence all the terms in any of the sums on the left hand 
side of equations (2) or (3) are of the same order, i.e. the super-matrices A and P 
are partitioned in a conformable manner. 


Either of the sets (2) or (3) of m? simultaneous equations is sufficient to 
determine the m? unknown sub-matrices P;,. Actually, for a fixed value of the 
subscript k there are, in the set (2), m matrix equations, which are sufficient to 
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determine the m sub-matrices P;, (i=1, 2, ..., m) which constitute the k" super- 
column of P, independently of the remaining super-columns. An actual solution 
may be achieved by a process of elimination similar to that used for sets of scalar 
equations. It is, however, often more advantageous to use equations from both 
sets (2) and (3) at the same time. 


Another way of obtaining the sub-matrices of the inverse matrix would be to 
transform the given matrix into the product of an upper and a lower triangular 
super-matrix, both partitioned correspondingly, A=UL. Except for their diagonal 
blocks, these triangular super-matrices are uniquely determined. One of the super- 
diagonals may be chosen as a set of m arbitrary non-singular matrices of appro- 
priate order; the usual choice is a set of unit matrices. The next step would be to 
invert the two triangular super-matrices. The inverse of A is finally obtained by 
applying the rule for the inversion of a matrix product, A~'=L~'U~'. 


It is thus seen that some of the methods in use for the inversion of ordinary 
matrices may be applied, with appropriate changes, to partitioned matrices. 
Whether much is to be gained by partitioning is another question. When a digital 
computer of sufficient capacity is available, either elimination or iteration may seem 
to be preferable to partitioning, at least with a matrix which does not exhibit any 
particular structure. The advantages of partitioning become most apparent in 
special cases, and these call for special treatment. 


The analogies just described between operations with ordinary and with 
partitioned matrices should not cause us to overlook some notable differences. First 
of all, in the product of two matrices the order of multiplication is not indifferent, 
and therefore a scalar formula containing products cannot be directly translated 
into a matrix formula. Furthermore, while it is possible to form a determinant 
from the scalar terms of any square matrix, this operation has no direct analogue 
when rectangular sub-matrices are considered. There does not seem to be any 
possibility of giving a practical definition of an adjoint super-matrix (with a possible 
exception in the case of uniform partitioning), and so there is no analogue to the 
useful relation between the adjoint and the inverse of a square matrix. Finally, 
there is no analogue to Cramer’s rule. Transposition of the super-rows and super- 
columns will not in general yield a result identical with the ordinarily transposed 
matrix. 


Some remarks about the effect of symmetry of the given matrix with respect to 
its principal diagonal are in order here. As with ordinary matrices, symmetry will 
simplify and shorten calculations when partitioning is used. The sub-matrices of 
the principal diagonal are then symmetrical themselves, and the sub-matrices above 
the principal diagonal are equal to the transposed sub-matrices symmetrically placed 
below the diagonal: Aj,=A;,;._ The same holds for the sub-matrices of the inverse 
matrix, with the result that the number of sub-matrices to be calculated is reduced. 
These simplifications being quite obvious, we shall present methods applicable to 
non-symmetrical matrices and only mention the possible reduction in the number 
of operations if symmetry is present. In the complex field, similar remarks apply 
to Hermitian matrices. 


May 1957 161 


per- 
rder 
on- 
r in 
that 
ion- 
the 

(3) 
n of 
er a 
gt k, 
and 
9 the 
order 
hand 
nd P 
nt to 
f the 
nt to 


ERYK KOSKO 


4. Partitioning into 2x2 Sub-Matrices 


This is by far the most widely known and often the only practical way of 
partitioning a matrix which is to be inverted. It has been described by Duncan™’, 
Hotelling*’ and others and is discussed here as an introduction to the treatment 
of more complicated cases. 

Let the given matrix be 


A=\p BJ’ 


and assume its inverse to be partitioned in a corresponding way, 


P 


The four relations represented by equation (2) are explicitly 
BQ+CS=I (1x1, BR+CT=O0 (1 x2’), 
DQ+ES=O (2x1’), DR+ET=I (2x2’). 


From (1 x2’), transferring the second term to the right hand side and 
pre-multiplying both sides by B~’, gives 


R= —-B-'CT. ‘ ‘ ‘ (4) 


Substituting this into (2~2’), and collecting the pre-factors of T into one 
bracket, gives 


where F is an auxiliary matrix of the same order (n,) as E: — 


F=E-DB~-'C. ‘ ‘ (6) 
From equation (5) follows 
(7) 
which, substituted into (4), yields 
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Similarly, the sub-matrices which form the first super-column of P are found 
to be 


Q=G-', . : (9) 
S=-E-'DG-', . . (10) 


is an auxiliary square matrix of the same order (7,) as B. 


The procedure just described involves no less than four inversions, albeit 
performed on matrices of lower order than the given one. The number of inversions 
can be halved if relations of the set (3) are also used, as follows: — 


QB+RD=I 1), QC+RE=O (I’x2), 
SB+TD=O (2’~x 1), SC+TE=I (2’x22). 


Transferring the second term in (1’ x 1) to the right hand side, substituting for R 
its value from equation (8) and post-multiplying both sides by B~', gives 


Transferring the second term in (2’ x 1) to the right hand side, substituting for T its 
value from equation (7), and post-multiplying both sides by B~', gives 


S= —F-"DB-’ ‘ (13) 


The expressions (12) and (13) for Q and S, alternative to (9) and (10), are 
preferable because they do not contain any inverses other than B~' and F~', which 
were already used to obtain R and T. 


The whole computation may be arranged conveniently in the following way. 
First, compute the sub-matrices of an auxiliary super-matrix 


=| 
in which B= C=-B-C 
D=-DB-! E=F=E+Dc, cin 
or =E+DC 


This part of the calculation begins at the top of the principal super-diagonal and 
ends at its bottom; we call this a forward sequence of operations. Then, by a 
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backward sequence of operations we compute the sub-matrices of the final: matrix 


_fQR 
“i 
in which Q=B+CS; 


T=E-', 
The alternative expressions for E and for Q may be useful for checking purposes. 


This process involves the following operations: two inversions of sub-matrices 
(one of order n, and one of order n.), six matrix multiplications and two additions. 
With a symmetrical given matrix two of the multiplications are replaced by 
transpositions. 


It would be equally feasible to compute an auxiliary matrix by a backward 
sequence, using E=E-! and B=B—CE-'D (identical with the matrix G defined 
by equation (11) ); the final matrix would then be obtained by a forward sequence. 
The result is obviously identical in both cases, although the auxiliary matrices in 
general are different. These results are easily verified by substituting into the 
equations (i x k’) or (i’ x k) the values of Q, R, S and T thus obtained. 


The order in which the various sub-matrices appear in any of the matrix 
products occurring in the formulae is easy to remember by noting that it is the 
same as the order in which they appear in the original or auxiliary super-matrix 
when following a clockwise course. It will also be noted that the final matrix is 
obtained from the auxiliary one in a manner somewhat similar to that in which the 
auxiliary matrix is obtained from the original one. There is, however, a difference 
in the signs of the non-diagonal terms and in an addition taking the place of a 
subtraction in the last sub-matrix computed in each of the sequences. 


5. Inverse of an Augmented Matrix 


5.1 THE GENERAL CASE 


It is sometimes useful to present the inverse of a 2 x2 super-matrix as the sum 
of two matrices, P=P,;+P;;. These are 


CF-'D, [DI 
and | F-'D, = (19) 
with C, D and F defined by equations (6) and (15). This follows immediately from 
the expressions (7), (8), (12) and (13) for the sub-matrices of P. Both P, and Pn 
are. singular, the first being of rank n, at most, and the second of rank n, at most. 
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Thus, when the inverse of a matrix B is known, it is possible to obtain the 
inverse of an augmented matrix A in which B is a diagonal block, by adding a 
correction term equal to Pj;. 

5.2. INVERSE OF A BORDERED MATRIX 

This case occurs when one of the diagonal blocks of the augmented matrix is a 

zero matrix, E=O, say. In this case the matrix F defined by equation (6) becomes 


The sub-matrices composing P are then obtained as in the general case from 
equations (7), (8), (12) and (13), with the restriction that the order n, of the zero 
sub-matrix must not exceed that of B. Otherwise, the auxiliary matrix F is singular 
and the inversion is not possible; this is because, in the expression for F, the rank 
of F cannot exceed that of B-?. 


In the particular case of a uniformly bordered matrix (n,=n.), C and D are 
square matrices; assuming them to be non-singular, we have 


and Q=B-'-B-'=O. 


The inverse of the bordered matrix is then 


oO 
p=(0., « » @ 


5.3. APPLICATION TO PROBLEMS INVOLVING STATIONARY VALUES OF QUADRATIC 
ForMs. (METHOD OF LAGRANGE’S MULTIPLIERS). 


Let u=} . (22) 
i,k=1 
be a quadratic form of the n variables x;(i=1, 2, ..., m) which may be considered 


as the components of a vector x; and let the matrix B=[b;,] of the coefficients be 
symmetric and non-singular. Further, let the variables x; be connected by n’ linear 
relations (n’ < n) 


t=1 


Then the conditions for a stationary value of u are 
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where the A, are undetermined multipliers, considered as the components of a 
vector A. After performing the differentiations, these conditions, together with the 
relations ¢,=0, can be written, in matrix notation, 


le olla] 


where C is the x 7’ matrix [c,,] and C’ its transpose. 


In accordance with previous results (with changes in the signs of some of the 
terms), the complete solution of the problem reads 


x) _ [°] 
li ] [ F-"C’ -F-'JldJ’° 
where C=B-'C, C’=C’B-' and F=CB-'C. 
The value of the multipliers being of no interest in itself, only the result 


is retained. 


The stationary value of u then becomes 
uy =4x,'Bx, =}4 dF-'d. ; . (26) 


This type of solution is of interest in the strain energy analysis of redundant 
structures when the statically indeterminate quantities must satisfy supplementary 
conditions. 


6. 3 x 3 Partitioning 


The solution for the general case of 3 x3 partitioning can be derived by a 
repeated application of the procedure described in Section 4 for 2x2 partitioned 
matrices. The computation may be arranged in a scheme which involves two 
auxiliary super-matrices. The operations consist of: three inversions (of orders 
n,, n, and n,;), 24 matrix multiplications and 12 additions. In the symmetrical 
case, nine of the multiplications are replaced by as many transpositions. The detail 
of the calculations has been given in Refs. 14 and 15. 


The special case of uniform partitioning has been treated in Ref. 7. 


A case of interest is one in which the two sub-matrices at the extremities of the 
secondary diagonal are zero sub-matrices. The given matrix may then be written 


B C O 
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Assuming the inverse to be 
P=(P,.]. (i, k=1. 2, 3) 


the relations between sub-matrices corresponding to equations (2) and (3) are 
formed. The details of the calculation are left to the reader to work out. From 
equations (1 x 2’), (2x 2’) and (3 x 2’), the unknown sub-matrices P,, and P,, are 
eliminated, giving 


P,.=L-’, 
then 
an 
Finally, the sub-matrices not belonging to the principal diagonal are determined : — 
P= ~B P,, = — P,,DB™", 
P,,= -P,, P,.= 
P,, =B~'CP,, P,,= Ds. 


These computations fit easily into a two-stage scheme similar to that of 
Section 4. First, an auxiliary matrix is computed 


_ |B CO 
A=|D SF (29) 
O GH 
with B=B-}, C= -BC, 
D= DB, E=L=E+DC+FG F=-—FH, . (30) 
G=-HG, H=H-'; 
and the final matrix P is then derived with 
P,,=B+CP,,, P= 
P,y=P,.D, P,.=GP.., P,, =H+GP., 


Alternative expressions may be used for checking purposes: E=E+DC+FG, 


P,,=B+P,.D, P,;=P,.F. P,,—GP.,, and P,,=H+P,,F. 
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The computation consists of three inversions of sub-matrices, 14 multiplications 
and four additions. In case of symmetry, five multiplications are replaced by as 
many transpositions. 


This type of matrix arises when two main systems of linear relations 
(represented by B and H) are linked by means of a third system (represented by B), 
without direct coupling between the main systems. The method here described is 
particularly advantageous when the linking term is of low order, i.e. when n, is 
small compared with n, and n,. 


Reference 9 describes the solution of a system of simultaneous linear 
equations along similar lines. 


7. Inversion of a Matrix of the Continuant Type 


Consider a matrix of scalars having non-zero elements only in the principal 
diagonal and in the two sub-diagonals immediately adjacent. Every row or column 
contains only three non-zero elements, with the exception of the first and last row 
and column which contain only two non-zero elements. Such matrices are 
encountered in connection with systems of equations of the “ three-moment ” type. 
They arise in the formulation of ordinary linear equations of finite differences of 
the second order, and therefore are an important tool in the solution of a large 
class of physical and technical problems. The name “continuant” has been used 
for a particular form of these matrices, the determinants of which serve to form 
approximations of continued fractions. The evaluation of these is based on 
recurrence relations, and such relations will be established for the computation of 
the inverse matrix. 


Consider such a matrix of the fourth order 


a, b, 0 O 

b, 0 
& %% 

O° 0 


and assume the inverse matrix to be 
P=[pix} (i, k=1, 2, a 4). 
Multiplication of the i#* row of A by the k column of P will be denoted, as 


before, by (i x k’); the meaning of (i x k) is also clear, and the result in both cases 
is equal to 6,,. From (1 x 1’) and (2’ x 1) elimination of p., yields 
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Similarly, using the relations (3’ x1) and (3’ x2) to eliminate p,, and ps., 
together with (1 x2’), which permits elimination of p,., gives, from (2 x 2’), 


1 boc. 
Po2 >= 7 + 


bic, 


where f,=a,- 
a, 


In the same way, (4 x 1), (4’ x 2) and (4’ x 3) are used to eliminate p,,, p,. and 
P43; and (1 x 3’) and (2 x 3’) are used to eliminate p,, and p.,, so that (3 x 3’) is 
transformed into 


1 b,c; 
P33> 4 fs? 
where f,=a,— ‘ 


Finally, substitution into (4’ x 4) of the values of p,,, p,. and p,,, just calculated 
from (4’ x 1), (4 x 2) and (4 x 3), gives 


where 
3 


Now it is possible to obtain successively p,;, p2. and p,, from the recurrence 
relations already given. The remaining unknown elements p,, (i ~ k) are obtained 
from the equations which have served to eliminate them. Without going into 
details, the inverse of the continuant matrix can be written 


(33) 
CiCo83 8s _ bs 


bcs , 


ffs” 


The whole computation for a matrix of order n may be summarised as follows. 


1+ =1+ 


b,c,g. 


where g,=1+ 
1J2 
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First, calculate the auxiliary quantities f;, starting with f,=a, and using the 
recurrence relation 


Next, calculate the auxiliary quantities g;, starting at the opposite end of the 
principal diagonal with g,=1 and using the recurrence relation 


Then the elements of the principal diagonal are 


The elements to the left of the principal diagonal in any row i are obtained from 
those lying immediately to their right by the relation 


=1+ 


The elements above the principal diagonal in any column k are obtained from those 
immediately underneath by the formula 


Similar expressions are derived in Ref. 10. 


8. Inversion of a Super-Matrix of the Continuant Type 


By analogy to the continuant matrix just discussed, we shall consider a super- 
matrix partitioned m xm with at most three non-zero sub-matrices in each super- 
row or super-column, with the exception of the first and last super-row or 
super-column (which contain only two non-zero sub-matrices each), the non-zero 
sub-matrices lying all in the principal super-diagonal and in the two adjacent 
diagonals. As before, the blocks in the principal diagonal are all assumed to be 
square, but not necessarily all non-singular. Let the given matrix be written 


O 
D, B, C, O 
O D, B CG, 


in 


0 0 ... Rey Be Gi 


and let its inverse be denoted by P=[P,,} (i, kK=1,2,...,m) 
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The derivation of the formulae in this section is similar to that in the preceding 
one, and the details of the calculation are not repeated here. The following 
auxiliary matrices are first computed : — 


:=B,, —-DB,-'C,, 


(40) 
C,= B;-'C,, D;= (i=1,2 


All the B; matrices are assumed to be non-singular. It is convenient to arrange 
these matrices into an auxiliary super-matrix 


A= 


O 


which matches A in the arrangement and order of the sub-matrices. 


The final (inverse) super-matrix is then computed, beginning with the bottom 
diagonal block 


The remaining sub-matrices are then obtained by recurrence, proceeding 
backwards. In the principal super-diagonal 


To the left of the principal diagonal 


Pie (i > k). 
Above the principal diagonal . 


The expressions (42) to (45) represent a complete solution of the problem. The 
total number of operations is: m inversions of sub-matrices, of order respectively 
N,, Ng, ...5 Nm; 3(m—1) matrix multiplications necessary to obtain the auxiliary 
super-matrix, and (m?— 1) multiplications for the final matrix, i.e. a total of (m — 1) x 
(m+4) multiplications; and 2 (m-—1) additions of matrices. With a symmetrical 
matrix A, the number of inversions and additions remains unchanged, but (m- 1) 
plus 4 m(m-—1), i.e. a total of 4(m—1) (m+2), multiplications is replaced by as 
many transpositions. 
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It can be seen that both the inversion by partitioning 2x2 described in 
Section 4, and that of the special 3 x 3 super-matrix of Section 6 are particular cases 
(m=2 and m=3, respectively) of the method described here. When substituting 
these values for m in the formulae of the preceding paragraph, the correct number 
of operations is found for these cases. 


The sequence of operations is to a large extent arbitrary. Any of the m! 
permutations of the diagonal blocks is admissible, as it can be obtained from the 
given super-matrix by simultaneous interchange of super-rows and super-columns; 
only in most of these arrangements, the non-diagonal blocks would not all lie in the 
immediate vicinity of the diagonal ones—that would require some adjustment in 
the scheme. The only essential restriction is that all of the auxiliary matrices B, 
must be non-singular. The example of the bordered matrix (Section 5.2) shows 
that it is possible to have one of them singular in one arrangement but not in 
another one. 


Matrices of the continuant type arise quite frequently in the theory of structures 
and in other applications of linear algebra. Most of their significant elements are 
grouped in and around the principal diagonal, while the elements away from that 
diagonal are zeros. Iterative methods of inversion have often been applied to these 
matrices, as they are by nature well-behaved and ensure quick convergence. It is 
by partitioning, however, that it becomes possible to take full advantage of the 
special structure of the matrix and of the presence of a large number of 
zero elements. 


9. Special Types of Symmetry. Use of Elementary Operations 
9.1. DouBLy SYMMETRICAL ARRANGEMENT 


First consider a matrix of order 2n=N with a doubly symmetrical arrangement 
of its sub-matrices, as follows 


BC 


Both B and C are square matrices of order n, not necessarily symmetrical them- 
selves. From considerations of symmetry it can be expected that the reciprocal of 


A has the same structure 
QR 


To determine the sub-matrices Q and R, the general method of Section 4 can be 
applied, using a suitable pair of equations from the sets (2) or (3). However, that 
would not allow us to take full advantage of the very special structure of the matrix. 


Instead, operations familiar from the elementary theory of determinants are used, 
namely addition and subtraction of multiples of rows or columns to other rows or 
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columns. In matrix language these are known as elementary operations and are 
equivalent to pre-multiplication or post-multiplication, respectively, of the given 
matrix by appropriate non-singular matrices, sometimes called elementary matrices. 
In the present case we shall use as multipliers elementary super-matrices derived 
from the unit matrix by inclusion of unit sub-matrices in suitable positions and with 
suitable scalar factors. 


A can thus be reduced to a diagonal super-matrix 


D=E’AE= 


(47) 


where the post-multiplying super-matrix is 


48) 


and the pre-multiplying super-matrix E’ is the transpose of E. The diagonal blocks 
of D are 
C=2(B-©. . . . 9) 


Inversion of the diagonal super-matrix is an easy matter 


On the other hand, both E and E’ being non-singular, there is a relation between 
this inverse and that of the original matrix 


which can also be written 


Thus the inverse of the given. matrix is obtained from the inverse of the diagonal 
super-matrix by transposing the operations which lead from the given matrix to the 
diagonal super-matrix. The explicit solution may be written 

Bee 


This confirms the anticipated symmetry in the arrangement of the sub-matrices. 


The whole computation is much simpler than the application of the general 
method; it consists of the inversions of two matrices of order n plus some additions 
and subtractions represented by the elementary operations of equations (46) 
and (51). 
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Note that this computation cannot be applied to a symmetric matrix 


unless C’=C, i.e. unless C itself is symmetric. 


9.2. CENTRO-SYMMETRIC ARRANGEMENT 


Next consider a matrix composed of 3x3 sub-matrices arranged in a 
symmetrical manner as follows 


BCD 
DCB 


Equal sub-matrices are found in diametrally opposed positions; this may be called a 
centro-symmetric arrangement, although the matrix A as such need not be centro- 
symmetric. Let B and D be square matrices of order n,, and F a square matrix of 
order n,, without making any assumptions as to their symmetry. 


I oO -I 
A multiplying matrix E,=|O 2I O (54) 
I oO I 


is used, where the I’s denote unit matrices of appropriate orders, and the O’s are 
blocks of zeros. Performing the elementary operation on A (i.e. pre-multiplying it 
by E,’ and post-multiplying it by E,) yields a reduced matrix 


D=E,’AE, = 


where B=2(B+D), 
C=4C, F=4F, 


The reduced matrix is somewhat easier to invert than A: — 


BC)" 0 
Oo O 


From this the inverse of the original matrix is obtained by applying the transposed 
elementary operation as indicated by equation (51). 
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9.3. CENTRO-SYMMETRIC MATRIX OF EVEN ORDER 


Let A, be a matrix of order 2m in which @=@on4,-i, 2n+:-x- By partitioning it 
in 2 x 2 sub-matrices, we may write 


where B and D are square matrices of order n, and the upsidedown letters € and d 
denote these same matrices with the order of their rows and columns reversed. 


A reversal of the order of the rows in a matrix B may be represented by 
pre-multiplication of B by a reversing matrix J in which all the elements of the 
secondary diagonal are equal to one, all the other elements being equal to zero. 
Similarly, a reversal of the order of the columns is represented by post-multiplying 
B by the reversing matrix. Finally, simultaneous pre- and post-multiplication by J 
yield a completely reversed matrix : — 


ad =JBJ, and conversely B=J@j. . ‘ 


Noting that J is orthogonal, and hence J~-'=J, we also have the identities 
BJ=Jd and AJ=JB, 


The aim is to reduce, by means of elementary operations, the centro-symmetric 
matrix A, to the form (46). This is achieved by applying a multiplying 
super-matrix 


the result of the operation being 


_fB DJ] [BC 
with C=DJ=Ja 


This can in turn be reduced to super-diagonal form, as in Section 9.1, but the same 
tesult is attained more directly by combining the two operations into one, i.e. by 
introducing A from equation (61) into equation (47): — 


The combined multiplying matrix is 


Equation (51) is again used to obtain A.~' from D-'. 
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9.4. CENTRO-SYMMETRIC MATRIX OF UNEVEN ORDER 


Let now A, be a matrix of order 2n+1 with Gj,=Qon4o-i, onso-1. It can be 
partitioned as follows :— 


¢ 


where B and D are square matrices of order n,& and C represent B and D with 
the order of their rows and columns reversed, f is a scalar element, C is a column 
vector, G is a row vector, and D and D9 are C and G with the order of their 
elements reversed. 


A, is reduced to the form (53) by using a multiplying matrix 
Oo O 
O J 


where I and J are both of the order n. Then, combining this with the operation 
indicated by equation (55), we have 


D=E’/AE, . . . . . . 66) 

with E,=K,E,=|O 1 2 2 . @ 
o I] jy o J 


The final inverse is obtained in the same way as before. 


9.5. INVERSION OF A CIRCULANT SUPER-MATRIX 


As a last example of special symmetry in the arrangement of the sub-matrices, 
consider a matrix partitioned as follows :— 


B C D 
A=|D B ‘ (68) 
C D B 


with B, C and D square matrices of order n. This is analogous to a circulant 
matrix of order 3, the determinant of which is calculated according to Aitken"? with 
the aid of the complex cube roots of unity. 
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Here a multiplying super-matrix 
I I 
E=|I «¢,I «,I ‘ ; . (69) 
I 


is used, where <, and <, are these cube roots: 


with the properties 


If E is used as a post-factor of A, the pre-factor shall be the inverse of E (or 
the conjugate of E, which is equal to 3E~'), rather than the transpose of the post- 
factor, as in the preceding examples. This operation reduces A to a diagonal 
super-matrix 


B OO 
D=E-'AE=|0 C (71) 
OOD 
where B=B+C+D, C=B+:,C+2,D, D=B+s,C+<D. (72) 


If A is a real matrix, the reduced matrix D is complex, and its inversion may 
present more difficulty than that of A. In some cases, however, especially when 
A itself is complex, it may be simpler to invert the three matrices B, C and D of 
order n than the given matrix of order 3n. We have then 


It would be easy to present further examples of symmetric arrangements of 
blocks within matrices which lend themselves to reduction to simpler form and 
thereby to easier inversion. 


10. Inversion of a Partially Modified Matrix 


When a matrix has already been inverted, it may be desirable to know the 
effect of modifying some of its elements on the inverse matrix. Assume that the 
modification is of an additive nature, i.e. that the modified matrix A is equal to the 
original matrix A, plus a corrective matrix A which contains a fair number of 
zero elements, 


A=A,+A. . ‘ ‘ ‘ ‘ . (74 
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Partitioning of these matrices permits separation of those areas affected by the 
modification from the rest. Assume then that the non-zero elements of A are 
grouped in the first m rows, the order of the matrices being n>m. Such a 
grouping can always be achieved by a suitable interchange of rows. 


Let A,~' be square partitioned into 2 x 2 sub-matrices 


Q being a square matrix of order m, and T a square matrix of order n—m. A 
corresponding partitioning of A is denoted by 


a=[8 (16) 


To invert the modified matrix A, it is advantageous to present it in the form 
of a product rather than of a sum:— 


Then the inverse is also obtained as a product 
Comparing equations (74) and (77), it is seen that the “ quotient ” matrix K is 
K=I1+AA,~—'. 


Putting for A and A,~' their partitioned values, from (75) and (76), gives 


(19) 


with M=I+BQ+CS, N=BR+CT. 


This is a triangular super-matrix, and its inverse is 


80) 


provided that M is non-singular, which we shall assume. Substituting this in 
equation (78) gives.the result 
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A simplification results from considering the difference between the inverse of 
the original matrix and that of the modified one. This correction may be developed 
into a matrix product 


or more explicitly 


The expression (82a) shows the correction term to be a singular matrix of rank 
m at most. 


Thus, instead of inverting the modified matrix of order n, it is sufficient to 
invert a matrix of lower order m, using the results of a previous inversion of the 
original matrix. 


Had the non-zero elements of A been grouped in the first m columns, 


a=[p ol 


a “quotient ” matrix would have been introduced, as in equation (77), but on the 
right hand side of A,, as follows: — 


A=ALL, 


with L=I+A,~'A. The resultant formula would read then 


with M=I+QB+RD, 


a square matrix of order m. 


These formulae become somewhat simpler when only one of the sub-matrices 
of A is significant. The case B+ O,C=D=O, has been discussed by the author”) 
in connection with structural stiffness matrices. The solution for the case of a single 
row or a single element being affected by the modification has been known for a 
long time. A recent contribution by Brock''*) presents the argument in very 
practical terms. 


Further cases of modified matrices are those of the augmented and bordered 
matrices, discussed in Section 5. 
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11. Count of Operations. Comparison with Elimination Method 


According to von Neumann and Goldstine"* the inversion of a general type of 
square real matrix of order n by the method of elimination requires 4 (4n° + 3n? +n) 
arithmetical multiplications and 4 (n?+3n) divisions. If the matrix is symmetrical, 
the number of multiplications becomes 4 (n*+2n?+n), the number of divisions 
remaining unchanged. Some other operations as additions, subtractions and trans- 
positions are left uncounted. 


When counting the arithmetical operations involved when the inversion is 
performed by the method of partitioning, the same formulae are applied for the 
inversions of the partial matrices. For a matrix multiplication of a (p x qg) matrix 
by a (q xr) matrix, the number of arithmetical multiplications is pgr. When both 
matrices are square of order n, this number becomes n*. The numbers of partial 
inversions and matrix multiplications are those listed in the sections of this paper 


dealing with the particular type of matrix. 


in the paper. 


On this basis, Table I has been prepared for several types of matrices discussed 
The results are probably on the optimistic side, as partitioning would 
necessarily involve additional logical operations and a more elaborate programming 
than straight elimination. In Table I, n represents the order of a typical 


sub-matrix in a uniformly-partitioned matrix. 


TABLE I 


COMPARATIVE COUNT OF ARITHMETICAL OPERATIONS 


Matrix Method of Partitioning Method of Elimination 
Invers- s Number of arithmetical Number of arithmetical 
Number |__!0ns = operations operations 
Type | § | of 
= & | Multiplications Divisions Multiplications Divisions 
General 6 | 1 On? +3n? +n n2+3n +6n?-+n 
Symmetric | 2n | 2 X 2 3 4n°+2n?+n n2+3n 4n3 +4n2-+n +3n 
General 3n| 3X3 3 | n | 24 | 1:5n?+4-5n |54n3 + 13-5n? + 1:5n| 4:5n2+4:5n 
Symmetric | 3n | 3 X 3 | | 16°5n?+3n?+1-5n} 1-5n?+4-5n| 13-5n3+9n2 + 1-5n} 4:5n?+45n 
Continuant 
General 14 | 1-5n?+4-5n |54n3 + 13-5n2+ 1-5n| 4-5n?+4:5n 
Continuant 
Symmetric | 3n | 3 Xx 3 9 | +3n?+1-Sn| 1:5n?+4-5n| 13-5n3 + 9n2 + 
Continuant 
General |10n | 10X10] 10| m | 126 | 146n3+15n?+5n | 5n?+15n |2000n*+150n?+5n| 50n?+15n 
e.g. 100 | 10 x 10] 10 | 10 | 126 147,550 650 2,015,050 5150 
Continuant 
Symmetric |10n | 10x 10) 10 | n | 72 +10n?+5n | Sn?+15n | 500n*+100n2+5n| 50n?+15n 
e.g. 100 | 10x 10] 10 | 10} 72 78,050 650 510,050 _ 5150 
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The table shows that second-order uniform partitioning offers no sensible 
advantage over the elimination method in the symmetric case, but in the asymmetric 
case there may be some saving of operations when partitioning. (A 2 x2 partition- 
ing with n,=1 and n,=N-—1 results in approximately the same number of 
operations as direct elimination. This is understandable, as in this case the two 
procedures are closely related). A similar situation, but with more pronounced 
differences, exists for third-order partitioning. 


With dontinuant matrices the advantage is always in favour of partitioning, 
rapidly increasing with the order of partitioning and much more pronounced when 
the matrix is asymmetric As an example, figures are given for a continuant 
matrix of order 100; they represent a saving of 84 per cent of the operations for a 
symmetric matrix, and 92-5 per cent for an asymmetric one. 


Special symmetries, as exemplified in Section 9, are obviously even more 
advantageous for partitioning, as the operations of reducing the matrix to super- 
diagonal form and vice versa are quite trivial. Only the inversions of the diagonal 
blocks contribute to the count of operations. 


12. The Diagonalisation of a Super-Matrix 


Summarising the findings of Section 9, the principal tool for inverting 
partitioned matrices is the elementary transformation of the given matrix into a 
diagonal super-matrix. The relation between these two matrices is an equivalence 
which, in the special cases considered, degenerates into a similarity or a congruence. 
The relation between the inverses of these matrices is the transpose of the relation 
between the originals. The problem of inversion is thus reduced to that of finding 
suitable multiplying matrices which will diagonalise the given matrix. In the 
examples discussed, these elementary transformations were not too difficult to 
construct. In the general case, however, diagonalisation is akin to the reduction 
of the matrix to a canonical form and may be far more difficult than the 
inversion itself. 


In the light of this process of diagonalisation how shall we interpret the 
methods of Sections 4 to 8? Are the auxiliary matrices A related in some way to 
the diagonal matrices D of Section 9? 


To answer these questions, observe first that a 2 x 2 partitioned matrix 


BC 


A=Lp 


can be reduced to upper-triangular form through pre-multiplication by a lower- 
triangular super-matrix 
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A similar reduction of the upper-triangular super-matrix 


lo 
OE 
is achieved through post-multiplication by another upper-triangular super-matrix, 
O EJLO I ~ tO EJ° 


Combining these two transformations leads to the diagonalisation of the 
given matrix:— 


a relation which can be abbreviated to 


where A is the given matrix, D the reduced diagonal super-matrix and L and U are 
lower- and upper-triangular super-matrices, respectively. The transposed relation 
between the inverses is 


P=A~'=UD"'L. ‘ . (84) 


Upon substitution of the values and after performing the multiplications, this agrees 
with the results represented by equations (16) and (17). 


It is seen now that the auxiliary super-matrix A defined by equations (14) and 
(15) is not in itself the result of transforming the given matrix in the proper sense 
(i.e. the result of some legitimate operations such as multiplications by suitable 
matrix factors). It is rather a compilation, a superposition of various blocks which 
are used at this stage of the calculation, namely B-', which is a part of D~', of 
E=E-DB-~’C, which is a part of D, and of the auxiliary quantities C=-B-'C 
and D= —DB~ which appear in the factors L and U. The auxiliary matrix thus 
combines in their proper place all the non-trivial quantities involved in this phase 
of operations; when properly interpreted, the auxiliary matrix presents a graphical 
plan of operations. 


A similar interpretation can be attached to the auxiliary super-matrix of 
Section 6. There, the process of diagonalisation of a 3 x3 partitioned continuant 
super-matrix is a little more complicated, but it also can be represented by multi- 
plications by appropriate factors 


MAN =D, 
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I O O 5 
with N=/|0 I O}. 
O O I O -H'G I 


Here again, the auxiliary matrix A, defined by equations (29) and (30), consists 
of parts of D, of D~’, and of the essential parts of the multipliers M and N. 
The final matrix is obtained as P=A-'=ND~'M, in agreement with the 
expressions (31). 


More generally, if we wish to eliminate a non-diagonal block A;, from a matrix 
A=[Aix], i.e. if we wish to transform that matrix in such a way that in the trans- 
formed matrix A a block of zeros appears at the place of A, this can be achieved in 
the following way. Pre-multiply the given matrix by a unit matrix to which a block 
— AixAj~* is added in the i super-row and k™ super-column. The transformed 
matrix differs from the original one in the blocks of the i super-row, a typical 
such block being 


Au = Ay — *An. 


If we now wish to eliminate from A another block in the same it" super-row, 
the result will generally be that non-zero elements appear again in the block 
previously eliminated. This explains why the method described in Section 6 is 
unsuccessful in dealing with super-matrices other than of the continuant type. 


We are now in a position to interpret the operations of Section 8 in terms of 
elementary transformations. Each step in the chain of “forward” recursions is 
equivalent to multiplication of the preceding result by a lower triangular pre-factor 
and an upper triangular post-factor 


Ais = L,A\U;. 
This begins with the given matrix A=A,, the first step being 


A,=L,A,U,, 


Each step eliminates two non-diagonal blocks C; and D; from the matrix. The last 
step leads to the diagonal super-matrix A,,=D which can be written 


D=LAU, 


with L=Ln-_,...L,L,, U=U,U, ... U,-;. 
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The “ backward ” process of obtaining the final matrix from the inverse of the 


diagonal one corresponds to the transposed continued multiplication 


... 


For practical calculations it is preferable to follow the step-by-step procedure 


as outlined in Section 8. The auxiliary matrix of equation (41) serves then as a 
graphical guide for the successive operations. 
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A Solution of the Flutter Determinant on a 
General Purpose Electronic Digital Computer 


D. B. GILLIES 

(National Research Development Corporation) 
and P. M. HUNT 

(The de Havilland Aircraft Company Ltd.)* 


Summary: A method is given for the solution of the flutter determinant on a 
general purpose electronic digital computer. The method has been programmed 
for the N.R.D.C. 401 Computer and details of this programme are given. A 
quicker programme, applicable when only structural damping is present and 
there are two degrees of freedom, is also discussed. Appendix I discusses in 
detail the types of errors which occur when the method is applied and Appendix 
II considers the anticipated loss of accuracy which will be encountered when 
using the method in cases where the number of degrees of fredom is large. 


1. Statement of the Problem 


The differential equations which arise in flutter calculations are of the form 


d*q ( d \dq ( 
apt b+ Ve) at V2 q=0. (1) 


where a_ is a matrix of order (nx n) of inertia coefficients 
bis a matrix of order (n x n) of aerodynamic damping coefficients 
cis a matrix of order (n x n) of aerodynamic stiffness coefficients 
d_ is a matrix of order (n x n) of structural damping coefficients 
e isa matrix of order (n x n) of elastic stiffness coefficients 
q_ isa vector of order (n x 1) of generalised co-ordinates 
V, is the aircraft velocity. 


The problem is to solve equations (1) for a range of values of V, and in particu- 
lar to find the values of V, which make one solution of (1) have constant amplitude. 
These velocities, usually called the critical velocities, and denoted by V., will, in the 
case of the lowest critical condition, be the boundaries between velocities at which 
the solutions of (1) are unstable and those at which the solutions are stable. The 
frequency of this oscillation occurring at the critical velocitiy is also required. 
Further, it is desirable that any coefficient in the set of equations (1) can be easily 


*Now of Ferranti Ltd. 
Received December 1955. 
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changed and the solution of the new set of equations obtained. In this way the 
effect upon V,. of varying mass balance weights or circuit stiffnesses, and so on, on 
the aeroplane can be investigated. 

The differential equations (1) represent the result of n degrees of freedom of the 
structure being both mechanically and aerodynamically coupled together. It is also 
desirable to be able to solve similar sets of equations representing the result of less 
than n degrees of freedom of the structure being coupled together. These equations 
are obtained from (1) by omitting suitable rows and columns from the matrices 
a, b, c, d, e. The process which is known as selecting the modes of vibration is 
useful in investigations to find which modes are coupling to give flutter when once 
flutter has been established for a particular value of n. It is also desirable that the 
rates of growth and decay of the oscillations at velocities near the critical velocities 
should be easily obtained. 


Many electronic analogue computers have been designed to solve this problem, 
but the following method is one which has been programmed for a general purpose 
electronic digital computer. It has all the facilities for carrying out the solution to 
the problem outlined. 


NOTATION 
Main Text 
Dots over symbols refer to differentiation with respect to ¢,. 
n number of degrees of freedom 
a matrix of order (n x n) of inertia coefficients 


b matrix of order (n x n) of aerodynamic damping 
coefficients 


c matrix of order (n x n) of aerodynamic stiffness 
coefficients 


d_ matrix of order (nxn) of structural damping 
coefficients 


e matrix of order (nxn) of elastic stiffness 
coefficients 


qr a generalised co-ordinate 


q vector of order (mx1) of generalised co- 
ordinates 


V, — velocity of aircraft 
V,. reference velocity 

V_ velocity defined by V=V,/V, 
V. critical flutter velocity 

o relative density 

B 


imaginary part of the root with zero real part 
which occurs at critical velocity 


t time 
t, defined by t,=t/(V/V,) 
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matrices of order (n x n) defined by 
A=a; B=b; C=c; D=d/V,; E=e/V,? 


arbitrary constants 

roots of polynomial given by equation (5) 
variable of polynomial of equation (5) 
polynomial defined by equation (6) 


matrices of order (n x n) defined by 
X=A; Y=BV+D; Z=CV°+E 


arbitrary numbers 


fA) 

defined by g (A;i)= 

II (,-A,) 

j=0 

defined by 
+k, + + koa i= II (A—A,) 

j=0 


coefficients of polynomial of equation (8) 


8 


defined by /,, 


Routh Test Functions 
defined in equation (9) 
defined by q;=p;/py 


defined at the end of Section 2 


defined by /,,.=L,, 


real variable 

complex variable 

small positive numbers 

Chebychev polynomial of degree 2n 


coefficient in Chebychev 


polynomial 
degree 2n 


of 


maximum modulus of <7... (A) in the unit circle 


error in evaluating a determinant 


error polynomial having values E; when A=A, 


coefficients in error polynomial 


error polynomial caused by truncating J, ; and 


due to rounding 
187 
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coefficients in error polynomial R (A) 
approximation polynomial to f (A) 

defined in text just before Theorem 1 
defined by L;(A)=J,, + 
defined in the text just before Lemma 3 
defined by M;(A)=M,, +Man, j 
defined in the text just before Lemma 3 


defined in the proof of Theorem 1 
positive numbers defined by equation (5). 


2. Method of Solution 


For convenience equations (1) are transformed before being solved. A reference 
velocity V, is introduced and equation (1) written as 


+] b+ Jat +1 ct q=9, (2) 


If t,=t/(V./V,) and the resulting equation is multiplied by (V4/V,)*, we 
obtain 
so that, if A=a, B=b, D=d/V,, C=c, E=e/V,’, 


Va d’q_ .. dq 


then Ag +(BV+D) q+ (CV*+E)q=0, (3) 


which is the set of equations that are actually solved. 


The solutions of the set (3) are given by 


= re exp {(2,+i8,)t,}, (=l,....m. @ 
Here ¢,, are arbitrary constants and 2,+i8,, with s=1, ..., n, are the roots 

of the polynomial 
| AA? +(BV +D)A+(CV?+E)|=0 . (5) 


of degree 2n. 


If, for each value of V, the values of 2,;+i8, can be obtained. then the problem 
has been solved, since by inspection of the signs of 2, (for each s) the solutions can 
be seen to be either stable or unstable and at the critical velocity at least one , will 
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be zero. Further (z;+i8,)/V are the damping factors which will give the rate of 
growth or decay of the oscillation at velocity V. If, however, the coefficients of the 
polynomial of which «,;+i8, are the roots have been obtained and one is only inter- 
ested in the stability of the oscillation, Routh’s Test Functions may be formed and 
used to determine whether the roots have any negative real parts and so whether the 
oscillation is stable or not. If all of the Routh Test Functions are strictly positive 
then the oscillation is stable, but if at least one is negative the oscillation is unstable. 
The boundary between stability and instability is given by the penultimate Routh 
Test Function being zero. 


The method consequently consists in finding, for each value of the velocity V, 
the coefficients of the polynomial of which «,+i8, are the roots and in either forming 
the Routh Test Functions for the coefficients of the polynomial or finding the roots 
themselves. Since the routine for finding Routh Test Functions is quicker, these 
are usually calculated for a range of values of V until velocities near V.. are obtained, 
when the roots of the polynomial are evaluated. Thus we obtain the rate of growth 
or decay of the oscillation at speeds near V.. and the frequency of the oscillation 
at V.. 


To obtain the coefficients in the polynomial, the determinant can be expanded 
directly, but an indirect method of expansion is probably better since, in general, it 
will be quicker. A well known method consists in evaluating the determi- 
nant at 2n+1 points and fitting a polynomial of degree 2n to these values by using 
Lagrange’s Interpolation Formula. In this way, once the values of the determinant 
at 2n+1 points have been obtained, the coefficients in the polynomial follow by one 
matrix multiplication. 


To show this, 
let f(A)=| AA?+(BV+D)A+(CV?+EB)| . © 
and, with chosen V, let A=X, (BV + D)=Y, (CV? + E)=Z, 
sothat f(A)=| XA?+ YA+Z | ‘ (7) 


By Lagrange’s Interpolation Formula, for arbitrary A,, .. ., Asn, we have 


2n 
2n 
fH= 2% FO 
| IL @Q—A)) 
i 
2n 2n 2n 
isi ixi 
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For convenience, denote 2 (A), lon 
II (A;— Aj) 
j=0 
2n 2n 
so that fay= {Tt a-a)} ca, 
i=0 wh 
2n 
and, if II (A-Aj)=A*+K, ke + han is 
j=0 


J 
iti 
then f{(A)= = +k, + +Kan, i} 


2n 2n 
] +)2n-1 kK, [= Kon, 3 (Ad) 


con 
= + + + Don con 
where Po | = ko, 0 on (A,) 
So t 
Pan Kon, 0 Ken, an (Asn) 
Hence, if /, .= Res 
II (A,—A,) mat 
j=0 
iss 
then py,.... Pon are given by 
Po i= I, 2n 
Pi fA) 
(8) 
To find f (A,),.... f (Aen) from (7) consists in evaluating 2n+1 determinants of 
order ( x n) and when these have been found substitution in (8) then gives p,...., 
Pon. Since the elements of the / matrix are functions of A,,..., Asn only and A,,..., 
Ayn can be chosen once and for all, the elements of the / matrix need be computed 
once only. The elements of the / matrix differ for different n, however, so that 
several matrices have in fact to be computed, one for each n, and all have to be held Con 
in the store at any one time if the facility for the selection of modes is to be available had 
and to work quickly. | 
May | 
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Having obtained the coefficients p,, p,,..., Pon, Routh’s Test Functions r,,..., 
Ton can be calculated from the following set of recurrence formulae. 


If po... Pon are denoted now by p,'”,..., Don, then 


where p,” are given by 
| 
pY=p,,,9-™ 
and 9=0 if i=0 (mod 2) ’ (9) 
=1 if i=1 (mod 2) 
and =0. 


The method which has been programmed for finding the roots of the polynomial 
consists in repeatedly extracting quadratic factors from the polynomial. For 
convenience, 


let q= > q2>= Po qn= Pm > m=2n, 


Po dD Po 


so that the polynomial to be solved is 
A™ + A™ +4m=0. ‘ ‘ ‘ (10) 


Then, if A* + 4A + 8" is a factor of the left hand side of (10) and a first approxi- 
mation to this factor is A?+,A.+,, a better approximation is 


where 2, 
Cm-1Cm-3 
i= 
-2 Cm-1Cm-s 


and (i=1,2,...,m) 
(i=1,2,...,m—2) 
b,=c,=1, b.,=c.,=0. 


Consequently, lim and lim 8,=B® and the roots of A? + =0 can 
be found in the “usual way. 
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3. Details of the Programme 


The method described in Section 2 has been programmed for the N.R.D.C. 401 
Computer and the programme requires the values of the A, B, C, D, E matrices 
(suitably scaled) to be specified, and it is then controlled by the hand switches. Four 
hand switches are used to indicate which modes are to be selected by the programme, 
a further sixteen hand switches are used to specify the value of the velocity in binary 
coded decimal, one hand switch is used to determine whether the roots of the poly- 
nomial are to be calculated or Routh’s Test Functions evaluated and one further 
hand switch is used when it is desired to read-in new coefficients. 


The steps in the programme then are: — 
(a) Select modes according to the hand switches. 
(b) Print an indication of the modes selected and the values of V. 
(c) Find X, Y, Z and f(A,),...,f (Asn) from (7). 
(d) Find py,...., Pon from (8). 


(e) Calculate 
either the Routh Test Functions r,,...,/2. from (9) and print 
them, 


or the a,+i8, and print them. 


The programme can be repeated now with either a new value of V or with a 
different set of modes selected or with a new coefficient in any of the matrices 
A,...,E. A special read-in routine has been written which reads a new coefficient 
into the store and then prints an indication of which coefficient has been read-in and 
also prints its value. 


The matrices A, B, C, D, E are scaled initially by multiplying the r* column of 
each matrix by a suitable factor so that the maximum element obtained in all 
these columns is 0:01. 


This process is repeated for r=1, 2,..., and the same process is now applied 
to the s“ row of each matrix for s=1,...,. Such scaling does not alter the solu- 
tion of the problem, as can be seen from equations (1), and ensures that capacity is 
not exceeded during the evaluation of the determinants. The determinants them- 
selves are evaluated by the method of pivotal condensation, the binary point being 
fixed until the pivotal elements are to be multiplied together, when floating point 
working is used. 


For the particular set of A,,...,A.» chosen the / matrices for n=2, 3,4 were 
calculated on desk machines and converted so that if /,.=L,. 2's then 
| | <1/Qn+1). 


This enabled the matrix multiplication of [I,, .] by [f (A,)] to be carried out using 
semi-floating point arithmetic instead of true floating point, with a corresponding 
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TABLE I 
by by Difference ros 
desk machine N.R.D.C. 401 (per cent) 
each V 
n=4 (Quaternary) 1073 10758 0-26 1 min. 40 sec. 
n=3 (Ternary) 1290 1289-4 0:04 1 min. 3 sec. 
n=2 (Binary) 1256 1256°4 0:03 35 sec. 


saving of time. The calculation of the Routh Test Functions is also done by semi- 
floating point but, when the roots of the polynomial are required, the coefficients in 
the polynomial are normalised and true floating point arithmetic is used throughout 
the rest of the working. 


Solutions obtained using this programme have been compared with solutions 
obtained on desk machines and these comparisons are given in Table I, together with 
the time taken by the computer for each velocity V. 


4. Case n=2 with No Structural Damping 


The programme which has just been described can be used when n=2 and a 
structural damping matrix appears in equations (1). However, if structural damping 
is not to be considered and n=2, the critical speed can be obtained directly by the 
following method. 


If y=(V,,/V 4)’, the coefficients in the polynomial 


+ + pA? + + (11) 
are given by 
14, Ay 
Are 
‘ dos bs: Ars 
Coo bs, Co, Are Cae C21 Are 
Bi, + [Ci Bi, +] Ors ] 
Coz Co, Doe bo Cae C2, Doe 
r= Cis | | Cir ] Ci C12 ] 
Coo C21 Coo Co, 


If the coefficients in (11) are appropriate to the critical speed, then equation (11) 
has an imaginary root A=if, say. 
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Hence, substituting and equating real and imaginary parts, 


pre 
P.B* — + =0 Ho 
pre 
+ ps =0, 
: and elimination of 8? gives 
6. 
PoPs’ — PiP2Ps + =9, 
which, by (12), is a quadratic in y=(V,/V44)’. wit 


This may be solved to give the critical flutter speed. The frequency parameter 
and flutter frequency can be obtained from = / (p,/p,). 


The foregoing method has also been programmed for the N.R.D.C. 401 Com- 
puter. The computer requires the matrices A, B, C, E to be specified and produces 
the critical velocities, frequency parameters and flutter frequencies in 7 sec. Since 
this programme is so fast it would probably be necessary to change the elements of 
the matrices A, B, C, E automatically in the store of the computer. For this reason 
the present programme also prints the matrices A, B, C, E, from which the solution 
has been obtained. The additional printing time is 28 sec. 


5. The Solution of the Stability Equations for Free Motion with 
Fixed Controls as a Particular Case 


— 


It is interesting to note that the programme described for the solution of the 
set of differential equations 


Ag+(BV +D) q+(CV’?+E)q=0 


can also be used to solve the stability equations for free motion with fixed controls. 


In general, the equation of small deviation of the longitudinal symmetric type 
and of the lateral and anti-symmetric type can both be written in the form 


AG+Bq+Cq=Q. . . . . . 


In general QO, () 


Q. (t) 
Q; (1) 


are, for the longitudinal symmetric motion, respectively the longitudinal force, 
normal downward force and pitching moment produced by the operation of controls 
or some external agency and are functions of the time. 


For the free motion with fixed controls, however, Q@=0 and the equation (13) 
reduces to 


A'G+Bg+Cq=0. 
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In this case the coefficients do not vary with speed and so we can use the 
programme with V=0 and A=A’, D=B’, E=C’ for the solution of this problem. 
However, the roots of the resulting polynomial (i.e. the stability quartic) would 
probably be required for each solution. 


6. Facilities Available when Solving with a Digital Computer 


Apart from the increase in accuracy the main advantages in solving the problem 
with a digital computer are as follows. 


(i) 


(ii) 


(iii) 


(iv) 


(v) 


(vi) 


(vii) 


(viii) 


May 1957 


The damping factors at any velocity may be obtained and so the severity 
of the onset of the flutter investigated. 


The initial scaling of the matrices themselves can be done by the 
machine. 


Any new coefficients which are put into the machine can be printed 
when they are read-in. Further, the whole matrices A, B, C, D, E can 
be printed at any time. Consequently a complete permanent record of 
the solutions is obtained. 


Arithmetical checks can be built into the programme to ensure that the 
working is correct. For example, in the evaluation of the determinants 
the row-sum checking system can be used. 


An overall check can be obtained on the critical speed V, by evaluating 
the determinant 


| -AB?+CV.2+E+i[8 (BV.+D)] | 


of complex elements on the machine and ensuring that it is zero. 8 is 
the imaginary part of the root with zero real part which occurs at the 
critical velocity. 


Solutions to binary flutter problems without structural damping can be 
obtained extremely quickly (7 sec.), so that many cases can be 
investigated. 


Automatic height changes can be incorporated in the programme. This 
is achieved by causing the machine to multiply the b matrix by /o and 
obtain the critical velocity with these new matrices. This speed will then 
be the E.A.S. critical speed at the height appropriate to which the 
relative density « has been chosen. 


Automatic changes in V can be carried out by the machine. Thus the 
velocity can either be increased or decreased by suitable steps automatic- 
ally, or the machine can carry out an interpolation on fr2n»-, to find a 
better approximation to V, when r.n-, will be zero. 


195 


D. B. GILLIES AND P. M:. HUNT 


(ix) The A matrix will always be positive definite and a programme can be 
written which will check that this is so if it is suspected that errors may 
have occurred in the computation of the elements of the matrix itself. 


(x) The same programme may be used to solve the flutter determinant and 
the stability equations for free motion with fixed controls. 


> 
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APPENDIX I 
THE CHOICE OF A,, Aj, +5 Aon 
va 

Real A; have been chosen, since only half as much calculation is required for A; 
‘, real as for A; complex or pure imaginary: evaluating f(A) for A complex or pure Me 
12 imaginary requires four times as much calculation as for A real, and produces twice as wil 
is much information, namely f(A) and f(A) where A is the complex conjugate of A. tha 
if In determining p (A) by this method, the principal sources of error are :— Lei 
| (i) The determinants’ f(A,) are not computed exactly, but are subject to an - 
§ accumulated rounding error. That is, f(A;)+¢,; is obtained, where ¢, is the Pre 
iz error in evaluating the determinant. Sin 
2 (ii) The round-off error in computing each scalar product p;= 3 I, jf (A;). This He 

error is small. 

Th 
(iii) Loss in the number of significant figures during the summation p;= & J, jf (A;) The 
: due to the subtraction of nearly equal quantities. This occurs when the ide! 
interval A,—A,, is too small or too large for the polynomial in question, so fa 
some terms p,A?"-' make only small contributions to the sum p (A) at every poi 


point A; used, and such coefficients are not accurately determined. This 
potential source of error can be guarded against by correct scaling and by Pro 
knowing the approximate value of the zero of largest modulus. 


(iv) The polynomial obtained is p (A)+E (A), where E (A) is the error polynomial 


having values <, at the A;. If the A, are unwisely chosen, some combinations Sim 

of small errors ¢; produce an enormous error polynomial, and in any case we The 

should like strict bounds on this type of error. M = 

This section will deal with minimising errors of type (iv), and will show that under | z | 
certain assumptions the optimum choice of A, is tcos(iz/2n) if the A, are to be 

distributed over the interval -t<A;<r. For simplicity, suppose t=1. (The case of But 
different t will be discussed later, and t=4 is actually used.) For t=1, the optimum 
A, are the extrema in the interval -1<A< +1 of the Chebychev polynomial of degree 
2n, that is, the 2n- 1 local maxima and minima, and the two end-points +1 and -1. 
We require the following well-known properties of T,, (A), the Chebychev polynomial 

of degree 2n:— Thi: 
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(a) T.,(A) is of the form +. that is, the coefficients of 
all odd powers of A are zero, the coefficient of A?" is positive, and successive 
even coefficients have opposite signs. 


(b) The 2n+1 extrema in the interval (+1, — 1) are symmetrically situated about 
the origin, and include 0 and +1. The values of T,,, (A) at these points are 


(- 1)", (- 1+}, (- 1)"+?,...., (- 1)". 


The type (iv) error is an error polynomial E (A) fitted to the errors ¢; in calculating 
f(A). Suppose each ¢,; is less than or equal to some positive quantity ¢, in absolute 
value. For any one choice of X,, A,, ...., Asn, there corresponds a “ worst” error 
polynomial which we define as that polynomial with maximum modulus in the unit 
circle for | ¢,;| <¢, i=0, 1, ...., 2m. That is, for fixed A,, A,, ...., Agn, We consider 
that error polynomial E (A) which results from the most disadvantageous combination 
of possible errors, and use its maximum modulus as a measure of the efficiency of that 
choice of A,, ...., Asn. We shall choose that set A,, ...., A, for which this measure is 
smallest. Since p (A) is obtained to estimate the location of the (usually complex) zeros, 
and particularly whether they are all to the left of the imaginary axis or not, this measure 
of possible error is reasonable. In fact, the solution to this problem minimises an 
impressive list of other things as well. One that it does not minimise is the maximum 
value on the real line, but this is very small anyway. 


Method of Proof. +¢T., (A) are the worst polynomials for optimum choice of A,. It 
will be shown first that, for any other choice of A;, worse polynomials exist, and second 
that, for optimum A,, the worst polynomials are, in fact, +¢T7,, (A). 


Lemma 1. For any choice A,, A,, ...-, Aon, the polynomial ¢<T,,,(A) is among the 
possible error polynomials. 


Proof. Set ¢;=¢T 4, (A,). 


Since the extrema of | 7,,(A)| are +1, | T,.Q)|<1, so | =¢| Tan Q)|<e. 
Hence ¢; are possible errors. 


There exists a unique polynomial of degree 2n, taking values ¢,; at A,, i=0, 1, ...., 2n. 
The error polynomial defined by ¢, and ¢7,,, (A) both satisfy this definition, so they are © 
identical. 


Lemma 2. The maximum modulus M of ¢T7,, (A) in the unit circle is assumed at the 
points + /(- 1). 


=(-1)"(to +ten) 

The modulus of both numbers is (¢,+f,+....+ton), SO it is sufficient to prove that 
Or to prove that | for any 
|z|<1. 
But | e7,,, (z) | =e | + | 

SeE(t, | z +t, | + Fan) 

Since | z| <1. 
This completes the proof of the lemma, 
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Definition. The interpolating polynomial L,; (A) for the point A; is defined as 
Ww 
2n (rX-A, 
ki 
2n 2n 
So p()= =L,A) fA) and E (A)= (A) &. 
i=0 i=0 
Note that the 2n zeros of L, (A) occur at the points A,, Ki, so L;(A)0 at all other 
points. 
of Theorem I. If not all A, are situated at extrema of T,,, (A), then there exists an error 
; polynomial with maximum modulus greater than M. D 
Proof. To prove the theorem, it is sufficient to construct an error polynomial whose Le 
lie modulus at ,/(-— 1) is greater than M. By Lemma 1, ¢7,, (A) is one error polynomial. Pr 
Since not all A; are situated at extrema of T,, (A), there exists some i (called i,) for pr 
which ¢ < yn <E. as 
i Hence the inequalities -¢<=8+¢T,, (Ai,) <e hold for at least one positive 5 and at Le 
: least one negative 6. Define an error polynomial E’ (A) by 
=eT,, A)=e; if ji, 
= (A,,) + é. an 
Then E’ (A)= eT (A) + 8Li, (A). Pr 
is Define a, b by a+/(-1)b=L;, {¥(-1)} and choose 5 positive if eT,, {/(-1)} 
is has the same sign as a, and otherwise choose 6 negative. Since | «T,, {/(-1)} | =M, If 
it follows that 
and, since a+ 1) b0, a?+b?>0 If i 
Now | E’{/(- D}| =] eTan{ D} +8{a+ 1) | Bu 
It f 
= J { {/(- 1)} + bap since 1)} is real 
{ + { 1)} +8 (a +64) 
since | /(- 1)}|=M. 
By (1) and (2), the sum of the second and third terms in the surd is positive, so pos 
Since the modulus of E’{,/(- 1)} is greater than M, the maximum modulus of E’ (A) 
in the unit circle is greater than M. 
This completes the proof of the theorem. whi 
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We have shown that if the A; are not all situated at extrema of T,, (A), there exist 
worse error polynomials than ¢7,, (A). Now suppose that the A; are situated at the 
extrema of T,,, (A). 


Then 1=A, >A, >.... >A, =O>.... 


and Aj 


Call the product of the numerators M,(A) and call the product of the denominators 


Lemma 3. The sign of D, is the same as the sign of (- 1)’. 


Proof. Dj=[(Aj;-A,) Aja) which is the 
product of j negative factors followed by 2n-—j positive factors, so it has the same sign 
as (- 1). 


Lemma 4. If M;(A)=M,, A2"+ .... +Mon, j 


then Mons, j= AWM . . (3) 


and the sign of M,, ; is the same as the sign (- 1)*.... (4) 


2n 
Proof. M; (A)= II (A —A,). 
If iAj and i<n, the two factors (A-A,), (A-Agn_;) occur in this product. Since 


= 


(A Ay) (A = = (A A) A+A)=A? 
If i=j and jn, the factors (A A,), (A — Agn_;) occur. 
But (A Aj) (A— Agn_)=A—0) (A +A) =A? 
n-1 
It follows that M; (A)=(A? +AA,) IT (A?—A,*)_ if 
iti, 
n-1 
and M, (A)= II (A?-A,”). 
i=0 
A product of factors of the form (A? — 4,7) has positive leading coefficient, the coefficient 


of all odd powers of A are zero, and successive even coefficients have signs alternately 
positive and negative. 


Hence, if j4n, M; (A) has the form 


where qo, ++ +> are all positive numbers, 
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and, if j=n, 
M,, (A)=q,A" — + (— 


since A, =0. 


Comparing coefficients in either (5) or (6) with those in 


M; (A)=M,, Man, ; 


we have (- ox, 

and M,,, ;=(- which has the same sign as (-— 1). 

Lemma If L; A)=(,, + +lon, then . . (7) 
and /,, , has the same sign as (— ‘ (8) 


Proof. L; (A). 


1 


Equating coefficients of equal powers of 4X, 


1 1 
i 2 
By equation (3), M 5=AiMax, 
1 1 
Hence D, D; ox, 


j=Ajlox, j . . (7) 


The sign of /,,, ; is the sign of D; times the sign of M,, ;, which is (- 1) (— 1)*=(- 1)*+/ 


An error polynomial E (A) has the form 
E(A)=SL; (A) ¢;, where | ¢;|<e, 


To determine the worst such error polynomial, we first determine the worst 


Lemma 6. For |¢,|<e, i=0, 1, ...., 2m, the worst polynomial of the form 
Low, + ( where is that for which | | =« 
and ¢,=(- 1)‘c,. The maximum modulus is attained at ./(- 1), and has value f,, | ¢, |. 
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But the maximum modulus occurs on the boundary of the unit circle, where | A | =1, so 
the maximum modulus is 


max | lox, A+ | = MAX | loess, |) 


since both sums are real numbers, and the maximum modulus is therefore attained at 
either A= +l orA=-1. 


=| = Dox, | | = +Aj) |. 


Since | A; | < 1, 1+A, is always greater than or equal to zero. The maximum sum is 
attained when all <,/,,, ; have the same sign throughout the summation. But, by Lemma 
5, the sign of /,, ; is (— 1)*+/ which is (— 1) times the sign of J, ;, so the sign of <; should 
be (— 1)/ times the sign of ¢,. 


Since the terms in the sum are then all of the same sign, the sum is largest when | ¢; | 
are as large as possible, that is, when | <;|=<«, i=0, 1, ...., 2n. 


Hence | «,| =< and <¢,=(- 1)‘e, for maximum modulus. But since A, are all situated 
at Chebychev extrema, <, (- 1)" T,, (A) =«, (— 1)'=¢, so the error polynomial in this 
case is ¢, (- 1)" T,, (A). 


Since T,,, (A) is of the form t,A7" +.... where t,, 
are all positive numbers, comparing the polynomials gives 


The maximum modulus of this polynomial is ¢,, | ¢, | and the value of the polynomial 
atA= 1) is (- 1)" + toy le, |: 


This proves the lemma. 


Lemma 7. For | ¢;|<«, i=0, 1, ...., 2, the largest positive or negative values of 
lan, Occur when | and ¢,=(- 1)‘e,. The value is then 


The proof of this lemma is similar to the proof of Lemma 6, but simpler, and 
will be omitted. : 


Theorem 2. If the A; are situated at the extrema of T,,, (A), the worst polynomials are 
+éT,, (A), and the maximum modulus on the unit circle is | ¢T,, { /(- 1)} |. 


Proof. The error polynomial E (A) is the sum of the n polynomials referred to in 
Lemma 6, and the constant coefficient referred to in Lemma 7. But the worst values 
of these polynomials occur for the same values of ¢,;, and the functional values at ./(— 1) 
are then all real and have the same sign as (- 1)" <,. Hence the maximum modulus of 
the sum is the sum of the individual maximum moduli and occurs at ./(- 1) when 
|e, | =e and ¢,=(-1)'¢,. The error polynomial is then ¢,7,, (A) which is +¢T7,, (A) 
and the maximum modulus is | ¢T,, {/(- 1)} |- 


Remarks. 
(i) Combining Theorems 1 and 2 it is seen that the unique optimum choice of A, is at 
the Chebychev extrema. 


(ii) This choice of A; has been shown to minimise the maximum modulus of the error 
polynomial in the unit circle. It can be shown that it also minimises the maximum 
value for any error polynomial of 


(a) the modulus of any coefficient of an even power of A, 
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(b) the sum of the absolute values of the coefficients, 
(c) the modulus of the functional value in any circle | R | <1, 


(d) the modulus of the functional value at z on the imaginary axis | z|<1 or 


(iii) The problem of finding optimum A, for fitting a polynomial of odd degree 2n+1 
has the same solution, namely A; should be at the extrema of T,,,, (A). 


(iv) With these results it is easy to establish several apparently new properties of 
Chebychev polynomials: Of the real polynomials less than or equal to 1 in absolute 
value in the interval — 1 to 1, the Chebychev polynomial has 


(a) maximum absolute value in any circle or at any point on the imaginary axis, 


(b) even/odd coefficients of maximum absolute value for polynomials of even/odd 
degree, 


(c) maximum sum of absolute values of the coefficients. 


(v) If the interval is —t to ¢ instead of —1 to 1, the entire argument can be repeated, 
substituting always A/t for A and for 


(vi) The value of | T,,{/(- 1)}| is $(1 + /2)?"+4 (1 - 2)". For large n this is about 
4(1+4 /2)*" or 4(2-414)?". The loss ot accuracy, measured in binary digits, is about 
(2:44n— 1) which, for n=4 is about 9 binary digits. This represents the largest error 
in estimating a functional value near a zero, supposing that that zero can be as far 
away from the origin as t\/(—1). The error is considerably less if the zero of largest 
modulus lies on a smaller circle, for example if all zeros lie inside a circle of radius 4, 
the loss in accuracy is only about (1-4n— 1) binary digits. In other words, the range 
—4=<) =<}, which is used after scaling, corresponds to some real range of damping 
factors. From the point of view of error polynomials, this range should be as large as 
possible, the limitations being that too large ranges do not determine sufficiently 
accurately the coefficients of low powers of A in p (A) because of type (iii) error. 


(vii) Analogue computation and this method suffer theoretically from the disadvantage 
that because of errors it is not proved rigorously that all roots lie in the left-hand 
half-plane. The following method remedies this situation. p(A)=f(A)+E(A)+R (A), 
where R (A) is the known rounding error and error caused by truncation of the coefficient 
I, ; That is, we can assign small coefficients r; to a polynomial = r,A?"-' which are 
known to be a bound on this error. If E(A)=e,A""+e,A""1+e,A""-?+...., then it 
is known that | e,, |< eto, | Cox,, | Stax. If then, the Routh Test Functions are 
re-calculated, allowing for rounding errors and a variation of +(éf.,+ | ro, |) and 
+(Etox+ | oxy, |) for coefficients of the form p,, and p,,,, respectively, it can be 
proved rigorously that for all possible error polynomials E (A)+R (A) the result remains 
the same, if the Routh Test Functions are still satisfactory. 


(viii) If the maximum errors for the various determinants are known, and are not all 
equal, a sharper result is possible. For the choice of A; given previously, coefficients 
of even powers of A of maximum modulus are obtained by setting | ¢,| equal to its 
maximum value for all i and sgn ¢;=(— 1)‘ and pre-multiplying the vector ¢,,.. -» 
by the matrix [/, ;], and coefficients of odd powers of A of maximum modulus are 
obtained by setting | ¢, | equal to its maximum value, for all i, 


and sgne,=(—1) ifixn 
and sgn¢,=(—1)'+! if i>n 


and pre-multiplying the vector (¢,, ¢,,... gn) by the matrix 
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APPENDIX II 
EXTENSION TO HIGHER VALUES OF n 


This Appendix attempts to explain why it is believed that if the calculation is carried 
out in double precision it should be possible to go to n=12 and perhaps higher. (Since 
the method is mathematically sound it is evident that, given sufficient precision, we can 
extend n, but from considerations of machine speed it is preferable if no more than 
double precision is required). 


The present programme is limited by the capacity of the 401 Computer in the 
following ways :— 


(i) The 401 store is only large enough to give single precision calculation. 


(ii) Of the 32 binary digits in a word, about 8 are sacrificed in scaling the 
matrices so that no intermediate results during the calculations can overflow, 
and about 9 are lost due to the error polynomial. If errors in evaluating 
determinants are counted, it is clear that the coefficients of the polynomial 
will not necessarily have more than 3 or 4 significant decimal digits. The 
loss of about 20 binary digits of accuracy is inevitable with the method, and 
the number of binary digits lost changes fairly slowly as n increases. 


(iii) The coefficients p; may vary considerably in size. The transformation \=2)’ 
implies p)’=27"p, and Pyn’=Pon, SO it is important to choose the correct 
size of interval A, to A,,. This can be determined automatically by the 
computer, and once found will probably not need to be changed over the 
whole range of velocities. There was not room in the store for this, so the 
range had to be set in advance. If one has very few spare binary digits, the 
correct size of interval must be determined accurately, and may even vary 
from root to root. Thus a root very near the origin is determined mainly 
by the coefficients of low powers of A, while a remote root is determined 
mainly by the coefficients of high powers of A. 


Ill-conditioning can always be detected by carrying out an estimate of the error 
during the calculation of Routh Test Functions, allowing p,, and p,,,, to vary by 
+t,,¢. If there is then uncertainty about the sign of one or more test functions, the 
offending coefficients {p;} can be found by repeating the calculation, varying only one 
of the coefficients (or all but one of the coefficients). 


The preceding discussion is intended to show that, no matter how pathological 
the polynomial is, it is possible to sense this fact, transform it into a non-pathological 
polynomial, and prove that the roots have been correctly located. Thus a polynomial 
of degree 20, having coefficients of 1, A, A?, ...., A’® near 1 and a coefficient of A”° of 
about 10-?° could have a root whose modulus was about 10. However it is felt that, 
since the polynomial arises from a real physical problem, such ill-conditioning will not 
occur. Apparent ill-conditioning, which is actually due to the mishandling of polynomials 
of high degree can be avoided, it is thought, by the proper use of the theorems and 
recommendations in this paper. 
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Review 


“ The Structure of Turbulent Shear Flow”, A. A. Townsend, Cambridge University 
Press, 1956. 315 pp. Illustrated. £2. 


When, towards the end of his great treatise on Hydrodynamics, Sir Horace Lamb 
came to the problem of Turbulence he prefaced his account by the remark, “ it remains 
to call attention to the chief outstanding difficulty of our subject.” Discouraging though 
it may be, any book written today, of similar compass and stature to Lamb’s, could 
hardly do better than repeat that statement. 


What complicates the study of turbulence is the essential randomness of the motion 
and the impossibility of describing, even statistically, its initial conditions: because 
every turbulent motion originates from a laminar flow which suffers a sequence of 
instabilities triggered by disturbances of unspecified statistical character, and the scale 
of the initial flow is never completely obliterated from the subsequent motion. This is 
one property of turbulence which distinguishes it from other random physical 
phenomena, such as the molecular motions in a gas, of which successful studies have 
been made. Another property, even more marked in its contrast to molecular behaviour, 
is the interaction between the motions of different portions of the fluid—an interaction 
which extends over distances comparable with the length scale of the entire turbulent 
flow—and for this reason the kind of model which is so effective in introducing the 
notions underlying the Kinetic Theory of Gases, that of replacing the fluid by an 
assemblage of discrete lumps or billiard balls, is unhelpful and misleading in the case 
of turbulence. To be sure, one of the earliest attempts at a description of turbulence, 
the Prandtl mixing-length theory, borrowed this type of idea from Kinetic Theory, but 
it is based on hypotheses so demonstrably unrealistic that it is now quite discredited. 
Indeed, so complete is the repudiation of mixing-length arguments in Dr. Townsend’s 
book that the reader will not even find Prandtl’s name in the bibliography! 


The basis of the modern approach to turbulence was laid by Sir Geoffrey Taylor 
more than twenty years ago (a precursory paper in fact appeared in 1921) and the 
development of his ideas into what is now a tolerably self-consistent, although incomplete, 
theory owes itself in no small measure to the researches of his colleagues Dr. George 
Batchelor and Dr. Alan Townsend. Up till the present time much of their work was 
available in a large number of separate, and occasionally contradictory, papers, but a 
coherent account of the present status of the subject (which, of course, includes contri- 
butions from laboratories other than the Cavendish) has been provided by Dr. Batchelor’s 
monograph Homogenous Turbulence and, more recently, by Dr. Townsend’s monograph 
in the same series The Structure of Turbulent Shear Flow. The two books are related 
even more closely than by the bond of title, for Dr. Townsend assumes a knowledge of 
the behaviour of the small-scale components of turbulent motion, the small eddies, of 
which Dr. Batchelor’s book gives an excellent account. 


In The Structure of Turbulent Shear Flow Dr. Townsend is principally concerned 
with the larger eddies which contain most of the turbulent energy. He begins his 
account by describing how a turbulent flow can be analysed mathematically and is 
intrepid enough to give a physical picture of an eddy. Hitherto, the term “ eddy ” was 
cloaked in physical anonymity simply by identifying it with nothing more or less than 
a wave-number component of the motion, but Townsend’s model has the delightful 
advantage of being actually observable, such as in the flow of water past the pier of a 
bridge or in the many regions of separation from the banks of winding streams. He 
then considers the effects on the eddies of distortion, imposed either by the mean flow 
or by the action of larger eddies, and in this way is able to give significance to the 
concept of eddy viscosity (although it must be admitted that throughout the book he 
remains non-committal in his attitude towards eddy viscosity). The introductory chapters 
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REVIEW 


also contain careful descriptions of the ideas of “ similarity ” and “ self-preservation ” of 
turbulent flows which are later applied to flows with free boundaries, namely, wakes and 
jets. In later chapters turbulent flows bounded by solid walls are discussed and the 
characteristics of these flows with those of free turbulence are combined in the final 
chapters to describe the more complicated, hybrid, structure of the turbulent boundary 
layer. The latter discussion includes the effects of pressure gradient, surface roughness, 
and of flow in tnree dimensions, such as that near the edge of a plate of finite span, 
with its attendant pattern of secondary flow. 


In a subject as perplexing as turbulence with its interplay between theory and 
experiment—neither of which at any one time is altogether free from criticism and 
emendation—it is refreshing to have a connected and eminently successful account from 
a person so uniquely fitted to discuss both aspects: and the reviewer has no doubt that 
the book will find a more permanent place in the literature than as the “ target for 
criticism ” which is what the Preface so modestly claims for it. Its appeal to workers 
in the field of turbulence and allied subjects is assured and it will certainly provide a 
stimulus to and inspiration for further research. But this review appears in an Aero- 
nautical Journal and it would be misleading to suppress doubts concerning the suitability 
of the book as a guide to turbulence for Aeronautical Engineers. The average engineer 
will find it very difficult to read, not because of any obscurity on the part of the author, 
but owing to the inherent complexity of the subject and unfamiliarity with some of the 
analytical methods used. At least, if the going is hard, it is of some comfort to reflect that 
turbulence is one of the few branches of Fluid Mechanics in which we can all claim 
to be in company with Sir Horace! However, the reviewer is convinced that every 
Aeronautical Engineer faced with problems involving turbulence would be amply 
rewarded if he were to make a conscientious and determining effort to study this book.— 
P. R. OWEN. 


Errata 
February 1957 (Vol. VIII, Part 1) 
In the paper by A. W. Babister “ Stability and Response of Systems Satisfying a 


Second-Order Linear Differential Equation with Time-Dependent Coefficients” there 
is a misprint in equation (1) on page 78. It should read: 


d@ d 


The author has asked for the following alterations to be made on p. 80:— 
Equation (9) to read 


co 


The first terms on the extreme left under “ Coefficient of” to read 


“ 


In the last line but one “ Cua, to read “c,_,.. 
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Errata 


November 1956 (Vol. VII, Part 4) 


In the paper by H. T. Jessop, C. Snell and G. S. Holister “ Photoelastic Investigation 
on Plates with Single Interference-Fit Pins with Load Applied to Plate Only” the authors 
regret some arithmetical errors in Table IV on page 306. This should read: 


TABLE IV 
TENSION AT HOLE BOUNDARY ON VERTICAL SECTION 


Tensile stress at C Increments 
(units)* (units) 


‘3 


4 


Hole 
Diameter 


(in.) 


come 


Low High 
load 


Low 
load 


High 
load 


(i) Low interference 


- 0-10 
0-13 
0-53 
0-67 


— 0-27 
0:37 
0:30 
+0-40 


— 1-57 
1:17 
— 0-94 
— 0-83 


(ii) High interference 


0:40 
0-74 
0:18 
0:54 


— 1:34 
1-27 
0-64 
~ 1-40 


Low 
load 


0-24 
0-28 
0-19 
+0:20 


— 0-35 
0-55 
-0-11 
0-27 


The figures in italics are those which have been*changed. 


High 
load 


—0-69 
— 0-44 
0-29 
0-21 
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1:47 1-20 
1-30 0-93 
1-47 1:17 
1:50 1-90 

3-34 2:94 2-00 — 0°59 
3-27 2-53 2-00 0-48 
3-24 3-06 2-60 -— 0-20 
3-47 2-93 2-07 — 0-35 
3 
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